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ABSTRACT 


Atmospheric turbulence severely degrades images of astronomical objects. 
Providing images that accuratefy reflect the true nature of these objects is essential 
to their understanding. Several object recovery techniques exist within the field of 
speckle imaging that produce accurate representations of astronomical objects. This 
thesis provides an in-depth comparison of two such techniques, Knox-Thompson and 
triple-correlation. 

Through computer simulation, this thesis accurately compares the abilities of 

both recovery techniques to enhance turbulence degraded objects by exploiting the 

di&action-limited information contained in short exposure, or "speckle", images. The 

* 

simulation produced these images by creating an object and several phase screens 
which simulated the effects of turbulence. Together, the object and the appropriate 
quantity of phase screens yielded the required short exposure images. Application 
of the Knox-Thompson and triple-correlation techniques to identical sets of these 
degraded images produced the resulting reconstructed objects, their signal-to-noise 
ratios and their azimuthal RMS phase errors. Comparison of these three factors over 
several imaging criteria concluded that the superior phase recovery technique was 
triple-correlation. 
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I. INTRODUCTION 


Scientific research of astronomical imagery has been continual since the 
invention of the telescope at the beginning of the seventeenth century. Initially, the 
limiting factor on astronomical image resolution was the quality and size of the 
telescope optics. As technology progressed, larger telescopes improved to the point 
where atmospheric turbulence (hereafter, simpty turbulence) became the limiting 
factor on image resolution. 

Turbulence severely limits the resolution of long exposure images produced by 
ground-based telescope imaging systems. Removal of turbulence corruption to 
improve image resolution is an area of extensive research. Optimal telescope site 
location minimized the effect of turbulence, but did not produce the near diffiaction- 
limited images desired. Removal of turbulence distortion occurred through the 
development of statistical methods, called speckle imaging techniques, that produced 
near difffaction-limited resolution. The basis of these techniques was the assumption 
that turbulence remains essentially stationary during a short exposure image of the 
desired object. Though distorted, these short exposure images retain diffi-action- 
limited information. 

Determination of the object’s power and phase spectra are separate operations. 
Although the power spectrum can be directly averaged while retaining diffraction- 
limited information, the phase spectrum cannot. Instead, either the cross-spectrum 
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(the Knox-Thonpson method) or the bispectrum (the triple-correlation method) are 
averaged because difEraction-limited information is retained. Then the phase 
spectrum is recovered from either the cross-spectrum or the bispentrum. The 
combination of the power and phase spectra generates the object’s Fourier spectrum 
which, when inverse Fourier transformed, provides the recovered image. 

Several factors affect the quality of image reconstruction. These factors 
influence both phase recovery techniques and include the amount of turbulence, the 
size of the telescope, and the light level of the object. The randonmess of turbulence 
and the effect of random photon noise, make exact image reconstruction impossible, 
however, increasing the number of short-exptwure images in the averaging process 
recovers a better image. This thesis compares the Knox-Thompson and triple- 
correlation phase recovery techniques under several different imaging conditions to 
determine their ability to improve the signal-to-noise ratio (SNR) of a degraded 
object. 
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n. BACKGROUND 


A. RESOLUTION 

In the image reconstruction process, resolution determines image quality. 
Resolution is defined as: 

the process or capability of making distinguishable the individual parts of an 
object, closely adjacent optical images, or sources of light [Ref. 1]. 

From a Fourier optics perspective, resolution is proportional to the high spatial 

fi’equency content of the imaging system. The standard for resolution is based on the 

image of two point objects (binary star) viewed through a telescope. The image 

consists of two overlapping Airy patterns with intensity 

!(♦) ■ ' ( 2 . 1 ) 


where J^(r|r) is the first order Bessel function and t|r is the normalized spatial 
frequency. The separation of the two first order fiinges of the two Airy patterns 
constitutes the so-called Rayleigh limit of image resolution 


AS = ^ , 


( 2 . 2 ) 


(in radians) for a circular aperture where X is the wavelength of light and D is the 
telescope diameter [Ref. 2]. Telescope imperfections and turbulence prevent 
attainment of the theoretical limit for large apertures. Speckle imaging techniques 
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remove distortions due to these factors and produce images with resolution near the 
theoretical limit. 

B. TURBULENCE EFFECTS 

Turbulence produces temporal and spatial variations in atmospheric density, 
temperature, and ind«c of refraction. The turbulence in the atmosphere perturbs an 
image, such as an Airy pattern from a point object (star), producing an extended 
image referred to as a seeing disk. The perturbation randomizes the electromagnetic 
phase front of the object. This randomization produces angular spreading, image 
wandering about its centroid, and scintillation or "twinkling". Turbulence effectively 
reduces the telescope’s resolving power by randomly attenuating the high spatial 
frequencies. 

C HISTORICAL PERSPECTIVE 

Until roughly 1970, attempts at solving the turbulence distortion problem were 
limited to finding the ideal telescope site. Generally, the sites were high in elevation 
and at locations regarded as having long periods of atmospheric stability. Even with 
great care in site selection, the typical angular resolution obtained was approximately 
one arcsecond, the maximum resolution attainable with a 12 centimeter telescope. 
Although phase distortions from turbulence constrained resolution, construction of 
large diameter telescopes provided enhanced light gathering capability. 

In 1970, a technique developed by Labeyrie enabled recovery of near 
diffraction-limited image Fourier moduli [Ref. 3]. The concept that a long 
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exposure image was blurred by turbulence fluctuations due to phase spectrum 
blurring and not power spectrum blurring, provided the basis of Labeyrie’s technique. 
Labeyrie determined that a short exposure image, about 10 milliseconds, would freeze 
the disturbance yet still contain near diffraction-limited information of the object. 
Taking many short exposure images, calculating their power spectra, then averaging, 
enabled a diffraction-limited estimate of the object’s power spectrum to be made. 
Labeyrie’s technique allowed high resolution measurements of binary star separations. 
However, lack of object phase information prevented faithful image reconstruction. 

In 1974 Knox and Thompson developed a technique for retrieving the object’s 
phase spectrum [Ref. 4]. The method uses the Knox-Thompson (KT) 
algorithm to provide an estimate of the object’s phase spectrum using the same short 
exposure images required for the estimate of the object’s power spectrum. The KT 
method calculates the average cross-spectrum of the object in Fourier space to 
determine the object’s phase spectrum. Calculation of the cross-spectrum involves 
determining the average correlation of spatial frequency pairs displaced from each 
other by a small frequent^ differential. The average provides a statistical phase 
difference approximation from which the object’s phase spectrum is obtained. 

In 1983 Lohmann, Weigelt, and Wimitzer developed another technique for 
retrieving the object’s phase spectrum [Ref. 5]. This method, referred to as 
triple-correlation (TC), also uses short exposure images to estimate the object’s phase 
spectrum. The TC method calculates the average bispectrum of the object in Fourier 
space to determine the object’s phase spectrum. Calculation of the bispectrum 
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involves determining the average of a third order correlation that consists of a 
frequency point, a point shifted by an offset, and a difference term. As with the KT 
technique, the average provides a statistical phase difference approximation from 
which the object’s phase spectrum is obtained. 

Either the KT or the TC techm'que determines the object’s phase spectrum 
which is necessary to produce accurate recovered images. Combining the 
reconstructed power spectrum and phase spectrum produces the reconstructed 
Fourier spectrum, which when Fourier transformed, yields the recovered image. 
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in. THEORY 


A. TURBULENCE MODEL 

Both the KT and the TC techniques utilize short exposure images to recover 
the object’s phase spectrum. A method developed by Tyler and Fried of the Optical 
Sciences Company, simulates turbulence to produce the short exposure images 
required to test these recovery techniques [Ref. 6]. This method requires the 
following three assumptions: [Ref. 7] 

1. Turbulence is represented by a single phase screen in the pupil plane of the 
telescope. 

2. Turbulence is isoplanatic, that is, the distortion from turbulence and the 
imaging system is considered shift invariant over the entire image plane. 

3. The images are quasi-monochromatic. 

With these assumptions, a single short exposure image becomes a convolution in 
image space 

i(^) =©(.?) ♦ siJi) , ( 3 . 1 ) 

where i (^) is the short exposure image intensity, o (x) the object intensity, and 
s (Jl) is the instantaneous point spread function. The vector j? represents the two- 
dimensional orthogonal spatial coordinates x and y. Using the convolution theorem, 
equation (3.1), becomes a product in Fourier space 


7 






J(il) = 0(3) • 5(3) , 


(3.2) 


where J(3) is the Fourier transform of the short exposure image intensity, 0(3) 
is the Fourier transform of the object intensity, and 5(3) is the instantaneous 
incoherent transfer function. The vector 3 represents the two^imensional 
orthogonal spatial frequency coordinates u and v. 

The point spread function, s {51 ), and thereby the incoherent transfer function, 
5(3), represent distortions from both turbulence and the imaging system. The 
assumption of stationary turbulence is accmate for short exposure images. 
Consequently, an instantaneous distribution of random phases (phase screen) 
approximates the instantaneous distortion of an image by turbulence. An array of 
random numbers filtered by a power spectral density function and corrected for low 
spatial frequency under-representation can simulate this phase screen [Ref. 6]. 
Therefore, 

5(3) = P(XF3)e^*<‘”» , (3.3) 

represents the instantaneous coherent transfer function, where P(Af3) is the 
transfer function of the telescope, is the instantaneous turbulence phase 

screen, F is the focal length of the telescope, and X is the wavelength of light [Ref. 
7]. Finally, the auto-correlation of the coherent transfer function, 5(3), 

5(3) = 5(3) ★ 5(3) . (3.4) 

yields the instantaneous incoherent transfer function required for equation (3.2). 
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B. PHASE SCREEN PRODUCTION 

Three common techniques for producing turbulence phase screens exist. One 
technique, the Fast-Fourier-transform (FFT) method, creates an array of filtered 
white noise and inverse Fourier transforms the array to real space providing the 
phase screen. A second technique, referred to as the Karhunen-Loeve (KL) 
expansion method, uses the KL expansion with a basis of Zemike polynomials to 
represent the phase screen. The third, hybrid, technique referred to as the 
Karhunen-Loeve-Fast-Fourier-Transform (KLFFT) method, combines the best 
properties of both techniques and manufactures phase screens which most accurately 
represent turbulence distortions. 

1. Fast-Fourier-Transform Method 

The FFT method provides a rapid means of generating a phase screen. 

# • 

Initially, creation of a square array of Gaussian-distributed random numbers of unity 
variance provides a representation of the phase over the entire aperture of the 
imaging ^tem being evaluated. The array amplitudes are filtered in Fourier space 
radially outward from the origin by the square root of the Kolmogorov power spectral 
density function 

F„(a) = 0.1517ro*^® I u , (3.5) 

where | il | is the radial distance from the origin in frequency units, and is the 
coherence diameter [Ref 7]. The origin is set to zero removing the constant (DC) 
term from the phase screen before applying the inverse Fourier transform. Two 
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phase screens result since complex numbers in the array consist of real and imaginaiy 
parts which are entirely distinct and statistically independent. 

Though the FFT method generates phase screens rapidly, it has 
deficiencies. The FFT uses a finite number of discrete points, and consequently, high 
and low spatial firequency cutoff occurs. High spatial frequenqr cutoff is minor since 
most of the wave firont error induced by turbulence is of low spatial frequency. Low 
spatial frequent^ cutoff is more serious as it induces under-representation of low 
spatial frequencies producing wave front tilt, or centroid position errors. Therefore, 
the associated structure function of the FFT-produced phase screen does not 
completely represent the 5/3 power law turbulence structure function. 

2. Karhunen-Loeve Expansion Method 

The KL expansion method provides an accurate method of generating a 
phase screen. Random phases associated with turbulence can be expanded in terms 
of a series of orthogonal functions (t), 

. (3.6) 

The expansion coefficients r^, are uncorrelated Gaussian random variables which 
represent turbulence statistics. The orthogonal functions provide the proper spatial 
dependence, thereby allowing the random phase to be evaluated anywhere within the 
aperture. The above expansion is refened to as the KL expansion. For a finite value 
of j, the KL expansion is the optimum basis whose eigenvalues represent the energy 
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content of the expansion coefficients, (t), and the total energy is the sum of these 
eigenvalues. [Ref. 8]. 


Zemike expansion coefficients, for a random phase screen, are Gaussian 
random variables. Unfortunately, these expansion coefficients are correlated and 
cannot be used as a KL basis set directly. However, the Zemike covariance matrix 
(the matrix of expansion coefficients) is useful in determining the KL expansion for 
turbulence. 

Three properties justify the use of Zemike polynomials as a basis set for 
the KL expansion to determine wave front turbulence. Use of the Zemike 
covariance matrix provides the necessary random variables required for the phase 
screen. Additionally, each eigenvector of the Z^emike covariance matrix is the 
representation of the KL function in terms of the Zemike polynomials. Further, each 
eigenvalue of the Zemike covariance matrix is the variance associated with the 
corresponding KL expansion coefficient. Wave front error induced by turbulence is 
an outcome of the these properties. [Ref. 6] 

The eigenvectors and the corresponding eigenvalues of the normalized 
Zemike covariance matrix are found which obey the relation 

CSi * , (3.7) 

where C is the covariance matrix, is the i** eigenvector and is the 
corresponding normalized eigenvalue. The eigenvectors are normalized so each 
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element of the eigenvector indicates the amount of the corresponding Zemike 
polynomial that is contained in the i ^ KL function expressed as 

' (3.8) 

p 

where K^(p) is the i“* KL function, is the component of thei“ 
eigenvector, and Zp(p) is the p'^ Zemike polynomial. Hence, the random wave 
front error produced by turbulence, (i'), is 

♦ (f) = , (3.9) 

where Yi is a set of Gaussian-distributed random numbers, f is the distance from 
the origin, and D is the telescope diameter. 

Though the KL expansion method is accurate, deficiencies exist, and care 
must be taken in its use. Calculating wave front distortion using the KL method 
requires an enormously large number of Zemike pofynomials to achieve enough 
accuracy. The required number is proportional to {D/Zg)^. Additionally, numerical 
inaccuracies exist in evaluating Zemike polynomials of high radial order. Therefore, 
to achieve the accuracy desired, avoidance of Zemike polynomials of high radial 
order is necessary. 

3. Karhunen«Loeve>Fast>Fourier-Transform Method 

The KLFFT method combines the fast computational speed of the FFT 
method with the optimum low spatial frequency representation of the KL functions. 
This technique requires a phase screen to be developed by the FFT method. The 
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first five KL functions are produced and applied to this phase screen to compensate 
for the under-representation of low spatial frequencies. With this compensation, this 
combined technique closely represents atmospheric turbulence. 

The KLFFT method is quite powerful. Since only five KL functions are 
produced and applied, the total number of operations per phase screen is 
approximately double that of the FFT method alone. Therefore, this phase screen 
production technique is relatively fast. [Ref. 6] 

C IMAGE RECOVERY 

The image recovery techniques represent methods for providing the 
reconstructed image from a turbulence-distorted object by utilizing several short 
exposure images of the object and a nearby star. The Labeyrie technique recovers 
the object’s power spectrum [Ref. 3]. The power spectrum provides the modulus of 
the Fourier transform of the object, the first part of the complex quantity required. 
Both the KT and TC techniques recover the object’s phase spectrum. Inverse 
Fourier transforming the product of the modulus and the phase provides a 
reconstructed image of the original object. 

1. Power Spectrum Recovery 

Equation (3.2) represents the image intensity for a single short exposure 
image in Fourier space. The time average power spectrum of several short exposure 
images in Fourier space is 
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<ii(0)|2> =|o(0)|2 • <\s(a)f> . 


(3.10) 


where the term on the left is the time average power spectrum of the image, the first 
term on the right is the object’s power spectrum and the second term on the right is 
the incoherent transfer function. Several images of a star under similar imaging 
conditions as the object determine this transfer function, which is the time average 
of the instantaneous transfer functions. The object and the star need not be in the 
same isoplanatic patch as long as the second order statistics of the transfer function 
are the same for both sets of exposures [Ref. 9]. Solving for the object power 
spectrum, equation (3.10) becomes 


|0(u)|* 


^ <|J(g)|^> 
<|5(0)|2> 


which recovers the object’s Fourier modulus. 


(3.11) 


2. Phase Spectrum Recoveiy 

Both the KT and TC phase recoveiy techniques recover the object’s phase 
spectrum by using the cross-spectrum and bispectrum averaging processes 
respectively. To recover the object’s phase spectrum, the KT technique calculates the 
average cross-spectrum while the TC technique calculates the average bispectrum. 


a. Technique Analysis 

The KT technique is the simpler of the two phase recovery methods. 
Determining the cross-spectrum is at the heart this algorithm. The cross-spectrum 
is defined as the product of two quantities in Fourier space: an array point and the 


14 




complex conjugate of an array point which is shifted, in frequency, from the original 
array point by a small offset, Au. The cross-spectrum, €5(3), is 

€5(3) = X(3) ‘ X*(3*A3) , (3J2) 

where X(3) is the array point and X* (3* A3) is the complex conjugate of the array 
point shifted by the ofrset vector A 3. 

The TC technique is a more complicated form of phase recovery. The 
bispectrum is defined as the product of -hree quantities in Fourier space: an array 
point, the complex conjugate of an array point which is shifted, in frequency, from 
the original array point by an offset, and an array point which is a function of the 
offset only. The bispectrum, B5 (3), is 

B5(3) » X(3) • X*(3*A3) • X(A3) , (3.13) 


where X(A3) is the array point which is a function of offset only and the other 
terms are the same as defined for equation (3.12). 

From equation (3.2), the average cross-spectrum is 


<I(3) • r(3*A3)> = 

[0(3) • 0*(a + A3)l • <5(3) • 5*(3*A3)> . 


(3.14) 


where the first term on the right is the object’s cross-spectrum and the second term 
on the right is the average incoherent transfer function cross-spectrum. The average 
bispectrum is 


15 






(3.15) 


<i{m • j*(0+Aa) • j(Aa)> = 

[0(a) • o*(a+Aa) • ocas)] • 

<5(3) • 5*(0 + A3) • 5(A3)> , 

where the first term on the ri^t is the object’s bispectrum and the second term on 
the right is the average incoherent transfer function bispectrum. 

The average incoherent transfer function cross-spectrum and 
bispectrum contain distortions firom both the turbulence and the imaging ^tem. The 
average phases resulting from these calculations for turbulence are zero. 
Imperfections in the imaging system produce phases which are negligible for both the 
cross-spectrum [Ref. 10] and the bispectrum [Ref. 5] calculations. Therefore, 
the average cross-spectrum and bispectrum of the incoherent transfer function is 
assumed to be unity. Equations (3.14) and (3.15) reduce to simpler forms 

0(3) • 0*(a+A3) - <r(3) • r*(3+A3)> , (3.16) 

and 


0(3) • 0*(3 + A3) • 0(A3) « 

(3.17) 

<J(3) • J*(3+A3) • J(A3)> . 
b. Phasor Spectrum Recovery 

Direct recovery of the object’s phase spectrum is not possible since 
the cross-spectrum and bispectrum phases are only known modulo In. The recursion 
algorithm fails when the cross-spectrum or bispectrum phase estimates are not equal 
to their principle arguments if multiple estimates of a single object phase spectrum 
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point are used. To avoid this problem without losing information, the reconstructed 
object’s phasor spectrum is determined instead. The recovery process wfll henceforth 
be referred to as phasor reconstruction. 

An arbitrary complex number N in phasor notation is 

^ = I I , (3.18) 

where | | is the modulus of the complex number and is the phasor in which <)> 

represents the phase of the complex number. Solving for the phasor, equation (3.18) 
becomes 


e‘* • - 


N 

\N\ 


(3.19) 


where the subscript ph denotes phasor. Therefore, solving for the phasors of 
equations (3.16) and (3.17) by dividing each by their respective moduli gives 

• 0^(3+An) = , (3.20) 

and 

' Op^(AU) = i^(i 2 ,Aa) , ( 3 . 21 ) 

where 




<j(u) • r{a*Aa)> 
|<J(u) • r*(iI+A0)>| ' 


(3.22) 


and 


rSca.AO) . <J(a)-J-(g^-AO) , - _ j(Aa)> 

^ |<I(a) • J*(il+A0) • J(A0) >1 


(3.23) 
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aTic tenns J^(Q,AO) and 1^(0,AQ) are four-dimensional quantities, since 
several offset values may be used to determine estimates of the object’s cross¬ 
spectrum and bispectrum. Solving for the oSset-shifted object phasor and applying 
the complex conjugate operator to both sides, equations (3.20) and (3.21) become 


0^(a*Aa) 


' I^(a,AQ) ] 

) 


(3.24) 


and 


o^{a*Am 


( i^ia.Arn V 

^o^(3)-o^(Aa)J • 


(3.25) 


Equations (3.24) and (3.25) provide the phasor spectrum necessary for image 
reconstruction. However, use of these equations is limited to infinite photon count, 
short exposure images which are unrealistic and useful for computer simulations only. 

c. Photon Noise effects 

Compensating for photon noise allows image reconstruction firom low 
photon count images, though more of these images are required in the process to 
resolve the object. During the recovery process of low light-level objects, the 
introduction of photon bias occurs. This bias corresponds to the cross-spectrum or 
bispectrum averaging of photon events with themselves and contributes no useful 
information to the average. With the bias removed, equations (3.11) and (3.12) 
become 
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cs(a) = if(a) • j:*(u+Aa) - jr(AiJ) , 


( 3 . 26 ) 


and 


BS(a) = X(a) 'XUO + All) ‘ X(AS) + 2Np - 

|j:(a)i2 + iJf(0+Aa)i2 - \x{Aa)\^ , 


where -X*(AiD removes the cross-spectrum photon noise bias, 

-\X{0 + AQ)|^, and -|jr(Au)|^ remove the bispectrum photon noise bias, and Np 
is the photon count [Ref. 11]. Therefore, for realistic image reconstruction, 
equations (3.24) and (3.25) become 


0 ^(a+A 2 ) 


rXTBZAS 

■Lph 


{Q.Am 





and 


Opft(a+Aa) 


( i^^(a,Am y 


(3.28) 


(3.29) 


where 


jJCTBIAS /0 _ <J(u) • J*(2-*-Ag) - J*(Ag)> 

^ ' |<J(u) • J*(2+Au) - J*{Aij)>| ' (3.30) 


and 


-TCBZAS 

•Lpb 


(a. Am = 


<1(2) • J*(2+A2) - |J(2)|2- |J(2+A2)|*- |I(A2)|2+ 2N^> 
|<J(2) • r(2+A2) - |J(2)|2- 1 J (2 + A 2 ) 12 - |J(A2)|2+ 2Np> 


(3.31) 
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Equations (3.28) and (3.29), allow realistic image reconstruction. 


d. Phasor Spectrum We^iting 

Phasor spectrum weighting provides a means for enhancing desired 
phasor estimates, increasing reconstructed image quality. Both phasor recovery 
techniques benefit from phasor spectrum weighting. Weighting techniques suppress 
higher frequencies and enhance the reconstructed image’s SNR. The method 
presented is the weighted least squares estimation approach. Matson showed that 
this method obtained the best results of four approaches anatyzed [Ref. 12]. 
The method weights the object’s phasor spectrum with the SNR determined from the 
variance of the cross-spectrum or bispectrum as follows: 




a^dmlSSiQ)]) 


± (Re[SS,(0)»l - 


J1 - 1 


(3.32) 




n - 1 


(3.33) 


o * V®* ]) + o* (Jjn[55(u) ]) 


(3.34) 


SNR » ^ , 

o 

where SS(Q) represents either the cross-spectrum or the bispectrum. 


(3.35) 
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e. Recovery Process 

Determination of the object’s phasor spectrum is afforded by the 
recursive method. An estimate of the image cross-spectrum or bispectrum for each 
short exposure image spectrum element provides the necessary information to 
determine the object’s phasor spectrum. The estimates for each array element are 
determined recursively, from the origin radially outward to the diffraction limit of the 
imaging ^tem. This method capitalizes on the inherent property of the phasor 
spectrum SNR which decreases radially as a frmction of increasing spatial frequency. 
These estimates are averaged over all the short exposure images and are weighted 
utilizing equation (3.3S). This process results in the object’s phasor array determined 
by equations (3.28) or (3.29). Once the reconstructed phasor array is found in this 
manner, it is multiplied by the square root of the power spectrum provided by 
equation (3.11). Inverse Fourier transforming this spectrum to image space provides 
the final reconstructed image. 

D. PHASOR RECOVERY TECHNIQUE COMPARISON 

There are several distinctions between the two phasor recoveiy techniques, 
aside from the obvious differences of equations (3.28) and (3.29). Since the recursion 
technique utilizes all the values from unity to the offset value in the averaging 
process, the greater the offset, the better the reconstructed image. For the KT 
technique, the modulus of the offset vector, | dd |, is restricted to be less than the 
seeing limit, Tq/X, whereas for the TC technique, the offset is unlimited. This 
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restriction allows the TC technique to surpass the KT technique in image quality by 
using a greater offset that includes more phase information. [Ref. 13] 

Shift invariance becomes important for low photon count images. The KT 
technique is not shift invariant and requires shifting of the degraded short exposure 
images to the center of the image array by, 


and 


X = 



• Ij(x) 
Ij(x) 


(3.36) 


yj • Ij (y) 



(3.37) 


prior to the recovery process. Image centroiding minimizes linear phase ramping in 
Fourier space, however, centroiding accuracy decreases with lower light levels 
resulting from the randomness of the image. Consequently, since the KT technique 
is not shift invariant, it performs poorly for low photon count images because the 
centroiding accuracy is decreased. The TC technique is shift invariant, eliminating 
the centroiding requirement. Since the TC technique does not require centroiding, 
the inaccuracy of the centroiding process at low light levels does not enter into TC 
image reconstruction. [Ref. 11] 
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E. SUMMARY 


The short exposure image is the key to effective image recovery. The high 
spatial frequencies within the short exposure image contain diffraction-limited 
information, though the object’s phase spectrum is randomized. Recovery of the 
phasor spectrum provides the information required for true image reconstruction. 

Three techniques for image reconstruction provide the recovered image. The 
Labeyrie technique recovered the object’s power spectrum. Both the Knox- 
Thompson and the Triple-Conelation techniques recovered the object’s phasor 
spectrum. Each phasor recovery technique yielded an independent array of phasors 
that, when combined with their corresponding power spectrum and inverse Fourier 
transformed, produced reconstructed images for comparison. Simulation of 
turbulence-degraded short exposure images by the Karhunen-Loeve-Fast-Fourier- 
Transform method allowed computer simulated comparison. 
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IV. COMPUTER SIMULATION 


A. COMPUTER REQUIREMENTS 

Turbulence simulation and image recoveiy processing, by their nature, require 
enormous amoimts of calculations. The simulation presented in this thesis creates 
an object and a turbulence phase screen. The phase screen distorts the object 
producing a short exposure image. Several of these short oposure images are 
utilized in the image reconstruction process by using either the KT or TC recovery 
techniques. The phase screens and short exposure images are presented in the form 
of two-dimensional square arrays. The arrays produced are of dimension 64 x 64 
consisting of 4096 elements. Several operations are performed on each element in 
these arrays throughout the entire simulation process. As a result, the requirement 
arose for a fast computer to process the arrays and for a large random access 
memoiy (RAM) to store them during the process. 

A personal computer was used for the simulation process and the data 
reduction. The computer, a Compaq Deskpro 80386/20 with 16 megabytes of RAM 
and a Weitek 1167 coprocessor, provided ample speed and convenience as long as 
array sizes did not exceed 64 x 64. Standard Fortran 77 was the language used 
throughout the simulation. A Microway 32 bit NDP Fortran-386 compiler provided 
the speed, precision, and array processing capabilities required for the simulation. 
The construction of each short exposure image and its subsequent cross-spectrum or 
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bispectrum estimation, required approximately 28 to 34 seconds to process with this 
computer-compfler combination, depending on the imaging parameters (coherence 
length, photon count, short exposure image quantity). Surfer and Grapher software 
from Golden Software provided plots of image data, phasor spectrum SNR 
comparison data, and azimuthal RMS phase error (henceforth, simply phase error) 
comparison. 

B. SIMUIATION PROCESS 

The simulation process compared the KT and TC phasor recovery techniques. 
The necessity for two programs, one using the KT technique and the other the TC 
technique, arose from RAM limitations. To ensure accuracy, both programs were 
identical accept for the individual phasor recovery subroutines. Additionally, the 
imaging parameters and the random number seeds were identical for each 
comparison run to ensure production of the same phase screens and, hence, the same 
short ocposure images. With identical short exposure images and object power 
spectra, the only distinguishable difference between reconstructed images resulted 
from the utilization of different phasor recovery techniques. 

The simulation utilized the speckle imaging prcKedure. This procedure involved 
the use of several short exposures, from 25 to 1600, to remove the effects of 
turbulence by means of an averaging process. This process determined the object’s 
power and phasor spectra by calculating the autocorrelation and the cross-spectrum 
or bispectrum respectively, for each short exposure image. At the end, these values 
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were averaged and combined, providing the object’s Fourier spectrum. A separate 
program filtered and transformed this spectrum to image space, yielding the 
recovered image. The simulation process is shown in Figure 4.1. 

1. Object Production 

Construction of an object that provided adequate detail to test the 
resolution of the two phasor recovery techniques was essential. Three objects were 
designed to compare the two techniques over several different image parameters. 
These objects were created, scaled, transformed to Fourier space, and normalized. 

a. Object Creation 

The first step in the process created the object. The option to 
construct one of three objects was provided. The first object resembled a finite¬ 
dimensional astronomical body centered in the array. The body was a convolution 
of a Gaussian function and a circular pupil function, giving it the appearance of a 
smooth planet Since phasor spectrum SNR declines radially with spatial frequency, 
roimd objects provide a less than ideal choice for image recovery comparison. 
However, increasing the detail on the body provided a means to test the resolution 
capabilities of the two phasor recovery techniques. Seven Gaussian functions of 
various size and depths at random locations on the body provided craters. These 
craters gave the body the appearance of an asteroid. The randomness of the craters 
on the asteroid provided the additional detail to test resolution. 
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)EGIN SIMUIATIOt 


Figure 4.1 Simulation Process 
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The second object resembled a binary star. This object consisted of two 
delta functions placed symmetricaUy about the center. Initially,this object provided 
the ability to troubleshoot the program since the phase spectrum of an equal intensity 
binary star involves a square wave pattern that was easily recognizable. After the 
completion of troubleshooting, the binary star allowed a test of the ability of the 
phasor recovery techniques to resolve point objects at very low photon counts. One 
point had twice the intensity of the other to reduce any ambiguities brought about by 
symmetry. 

The third object resembled a star. The object consisted of a delta 
function at the center of the array. The short exposure image of a point source yields 
the instantaneous incoherent transfer function representing distortion from the 
turbulence and the imaging system. The transfer function allowed recovery of the 
object’s power spectrum. 

b. Object Scaling 

The object had to be scaled before use. Object scaling assured that 
the imaging parameters retained complete frequency information within the given 
array size and ensured that the size of the object was within acceptable limits. One 
such parameter was rQ, the coherence length, which was a measure of the amount of 
turbulence present in the atmosphere [Ref. 14]. For the simulation, values 
of 0.206 and 0.103 meters were chosen. Another parameter, |u|, was the offset 
value. This value represented the number of pixels (array elements), in Fourier 
space, contained within the coherence length. The offset maximized the averaging 
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of cross-spectrum and bispectrum while preventing loss of high frequency information. 
O&et values of two and four were chosen for the simulation. The corresponding 
width of each pixel in terms of coherence length was determined and divided into the 
telescope diameter to calculate the number of pixels retained by this telescope imder 
the conditions of the above parameters. The number of pixels is synonymous with 
the frequency cutoff of the diffraction-limited telescope incoherent transfer function 
(.D/X ). The image array size limits the frequency cutoff value to 





(4.1) 


where is the frequency cutoff and n is the array size. If the frequency cutoff is 
too large, frequency information is lost. 

After constraining the imaging parameters, the field of view and the object 
size were ascertained. The field of view (FOV) is equivalent to the reciprocal of the 
fundamental spatial frequency 


FOV = |i2( • 



(4.2) 


where (X/r^) is the seeing disk. In general, the allowance of one seeing disk width 
between the object and each side of the array provided for object distortion effects. 
The size of the asteroid was then maximized to provide the most visual detail while 
satisfying the above criteria. Verification that the binary star and the star obeyed the 
above criteria was sufficient for those objects. 
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c. Fourier Spectrum Determination 

Production of phase screens occurred in Fourier space, while object 
production transpired in image space. Proper image production required 
transformation to frequency space. A two-dimensional Fast-Fourier-transform (2D- 
FFT) algorithm transformed the object to Fourier space. The FFT provides a fast 
and accurate method of transforming a discrete frinction to Fourier space. The 2D- 
FFT employs a ID-FFT provided by Gonzalez and Wintz [Ref. 15]. This 
ID-FFT determines the discrete Fourier transform of a complex one-dimensional 
array of numbers. The 2D-FFT simply calls the ID-FFT for each row then each 
column of the two-dimensional object array. Testing of the 2D-FFr by transforming 
a normalized pupil function then inverse tramforming enabled comparison between 
the results and the ori^al function. With double precision complex numbers, the 
2D-FFT provided accuracy to ten significant figures. In addition to determining the 
object’s Fourier spectrum, this 2D-FFT provided the means for transformation, from 
image space to Fourier space and back, extensively throughout the simulation. 

d. Fourier Spectrum Normalization 

Normalization of the object’s Fourier spectrum produced the correct 
number of photons in the short exposure image. In reality, each short exposure 
image furnishes a photon count, however for the simulation, the some value was 
chosen for all short exposures used for each image reconstruction run. Dividing the 
object’s Fourier spectrum by the photon count normalized the spectrum. Since the 
phase map was normalized by its DC value, the photon count was the same for the 
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object and for the short exposure image. The normalization was impoitant for 
Poisson noise generation introduced later in the process. 


2. Turbulence Phase Screen Production 

Construction of a turbulence phase screen that correctly resembled true 
atmospheric turbulence allowed accurate phasor recovery technique comparison. 
Testing the KLFFT method of phase screen production showed it represented the 
actual 5/3 power law structure function closely [Ref. 16]. Therefore, this 
method was adapted to produce the phase screens in this simulation. 

a. Gaussian Random Number Array 

Each phase screen involved an array of random numbers. The 
random numbers represented the random phases produced by turbulence. A 
Gaussian distributed random number generator, provided by the subroutine Gauss, 
produced the required random numbers that represented the randomness of 
turbulence statistics. 

b. Filter Function 

The simulation used a filter function that represented the square root 
of the Kolmogorov power spectral density function, equation (3.5). This function 
represented turbulence statistics and filtered each array of Gaussian distributed 
random numbers to provide the turbulence structure function. The resulting array 
elements when put in terms of the Rytov approximation, modelled the 5/3 
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power law structure function, and hence provided a model for turbulence, though the 
low spatial frequencies were under-represented. 

c. Karhimen-Loeve Functions 

The simulation iised the first five KL functions. Each filtered array 
was inverse Fourier transformed and then the five KL functions were applied in 
image space to compensate for the low spatial frequency under-representation. This 
application involved more than simple multiplication. The inner product of the KL 

* 

functions and the filtered array provided scalars which expressed the amount of each 

individual KL function contained in the array. These amounts were subtracted from 

the array. The technique involved random numbers with variances equal to the 

eigenvalues associated with their respective KL function. Multiplying the random 

numbers by their corresponding KL functions and adding the result to the array gave 

« 

the corrected phase screen. The resulting array elements when put in terms 

of the Rytov approximation, accurately depicted the 5/3 power law structure function. 

An effect of the KLFFT method was the inclusion of tilt in the short 
exposure image. An option was given in the simulation which allowed the removal 
of tilt setting the first two KL functions to zero. This option allowed phasor 
recovery without the presence of tilt and was utilized as a criterion for comparison * 

of the phase recovery techniques. 
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d. Incoherent Tranter Function 

Use of the Rytov approximation determined the incoherent transfer 
function of the turbulence and imaging ^tem. The Rytov approximation, 

, (4.3) 

provides the coherent transfer function of the turbulence and imaging system, where A 
represents the amplitude (set to one) and represents the realization of the phase 
screen determined above. Calculating the autocorrelation of the coherent transfer 
function produced the incoherent transfer function. This calculation first required 
multipfying the phase screen array produced by the Rytov approximation by the pupil 
function of the telescope. Squaring the modulus of the Fourier transform of this 
array and Fourier transforming and normalizing the result yielded the incoherent 
transfer function of the true phase screen. 

3. Object Degradation 

The product of the phase screen and the object’s Fourier spectnun yielded 
the short exposure image spectrum. The inverse Fourier transform of this spectrum 
resulted in the required short exposure image. Since the object was normalized to 
the photon count of the short exposure image and the phase screen was normalized 
to unity, the final step in the degradation process was to apply photon noise to the 
short exposure image. The photon noise effect was added to the short exposure 
image by entering each image element into a Poisson distributed random function 
generator provided by the Poisson function in the simulation. The generator returned 
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a random number at each image point drawn from a Poisson distribution with mean 
equal to that value. As the photon count decreased from infinity, the randomness of 
the returned values increased which caused the grainy, photon noise effect With 
photon noise included, the short exposure image was complete. 

4. Centroiding 

The simulation supported centroiding the degraded short exposure in 
image space. If desired, the image array was shifted to its true centroid first by 
column, then row, based on the centroids determined from equations (3.32) and 
(3.33). Determining the centroid of the binary star object modified with equal 
intensity stars checked the accuracy of the centroid subroutine. The object initially 
had an arbitrary translation from its centroid. Centroiding this translated object using 
the subroutine then analyzing its phase spectrum ensured the subroutine performed 
correctly. Image centroiding worked well on high light-level images of the asteroid 
above 10^ photons. Images below this level had unevenly distributed intensities and 
centroiding was ineffective. 

5. Image Recovery 

The simulation separated image recovery into two distinct parts. First, the 
Labeyrie technique provided the object’s power spectnun. Use of this technique for 
both programs ensured uniformity in comparison. Second, the KT and TC techniques 
reconstructed the object’s phasor spectrum, each in separate programs. 


34 






a. Power Spectrum Recovery 

Recovery of the object’s power spectrum provided the modulus of the 
object’s Fourier spectrum. Estimates of the autocorrelation were calculated for each 
short exposure image of the object and the point source. Averages of these results 
yielded the average power spectrum. Normalizing both power spectrum arrays by 
their respective DC values provided an equivalent intensity basis. Dividing the 
object’s average power spectrum by the point source’s average power spectrum 
removed imaging system errors. The square root of this result multiplied by the 
photon count furnished the object’s modulus. 

b, Phasor Spectrum Recovery 

The object’s phasor spectrum, with its modulus, determined the 
object’s Fourier spectrum. Each short e}qx)sure image provides estimates of the 
either the cross-spectrum or the bispectrum, for various offset values. Several of 
these estimates, each with a different offset value, determined a specific phasor array 
element by averaging these estimates over all short exposure images, then over aU 
offset values. 

Recursive calculation, outward fi-om the origin of the image array, 
supplied the estimates, and hence the object’s phasor spectrum. The number of 
estimates which existed within an estimation circle of integer pixel radius equal to the 
whole part of the pixel distance of the desired point determined the maximum 
number of estimates and the offset value for each phasor estimated. The estimation 
circle began at the origin, and the estimates of the phasor spectrum points began one 
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radial pixel distance beyond that. The phasor at the origin had a value of one, since 
it had no imaginaiy part, and it was normalized. Each phasor spectrum point 
included more estimations as the recursive process proceeded, out to the maximum 
radial o&et value. The same estimates of the phasor spectrum were found for every 
short exposure image and then averaged. Averaging a set of estimates associated 
with a phasor spectrum point produced the corresponding phasor value for that point. 

Using the phasor spectrum of the uncorrupted binary star object 
modified with equal intensity stars verified the phasor recovery process. Comparing 
the phasor spectrum of the binary star before and after the recovery process using 
complex double precision numbers provided a match for all points to ten significant 
figures. 

The phasor recovery process required an extensive amount of 
calculation. The time for phase reconstruction was approximately 30 seconds for 
every short exposure image. The process was made less time consuming by invoking 
Hermitian symmetry. Hermitian symmetry dictates 

Oi.j = Oli , ( 4 . 4 ) 

therein allowing half the number of calculations to determine the full object’s phasor 
spectrum. 

At the end of the phasor recovery operation, the variance of the cross¬ 
spectrum or bispectrum estimates provided an SNR value for each phasor element 
from equations (3.34) through (3.37). The square of the SNR for each estimate was 
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then multiplied 1 ^ the corresponding estimate when the phasor spectrum points were 
calculated. This calculation provided a weighted least squares estimation of the 
object’s phasor spectrum. When combined with the recovered modulus, the object’s 
Fomier spectrum resulted. 

6 . Azimuthal Signal-To-Noise Ratio 

Averaging the SNR values of the cross-spectrum and bispectrum provided 
an average SNR value for each phasor element. This average SNR array was then 
averaged azimuthally, one radial pixel value at a time, from the origin out to the 
cutoff frequency to provide a SNR as a frinction of radial spatial frequency. This 
radial SNR provided one means to compare the two phasor recovery techniques as 
well as to determine the frequency at which noise overcame signal to enable proper 
filtering. 

7. Fourier Spectrum Filtering 

The final result of the simulation process was a weighted least squares 
estimate of the object’s Fourier spectrum. Before inverse Fourier transformation, the 
object’s Fourier spectrum required filtering. A simple rectangular low-pass filtering 
method, which truncated spatizil frequencies beyond the radial frequency where the 
azimuthal SNR was unity, determined the object’s Fourier spectrum. This filtering 
process was provided by a separate program which included a method to determine 
the phase error of the recovered image’s Fourier spectrum. 
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8 . Azimuthal RMS Phase Error 

Measuring the phase error of the two phasor recovery techniques 
produced another means for comparison. The object’s phasor spectrum determined 
its phase spectrum. Computing the phase spectrum modulo Iv, then subtracting the 
result from the object phase spectrum of the true object representation provided the 
array point phase error. The square of this error was determined then averaged 
azimuthally one radial pixel value at a time, from the origin out to the cutoff 
frequent^. The square root of this average provided the azimuthal RMS phase error. 

C SUMMARY 

The simulation process obtained the reconstructed image from several short 
oqx)sures images. The process p^odrced the desired object and the phase screens 
to make the short exposure images required for speckle imaging. With several short 
exposure images of the object, the Labeyrie technique recovered its power spectrum 
and the Knox-Thompson and Triple-Correlation techniques recovered its phasor 
spectrum. Low-pass filtering removed noise and the SNR and phase error 
calculations presented means for comparison of the two phasor recovery techniques. 
The inverse Fourier transform of the filtered spectrum yielded the recovered image 
presented in the form of a contour plot. 
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V. SIMUIATION RESULTS 


A. RECOVERY TECHNIQUE COMPARISON CRITERIA 

Seven criteria provided a basis for comparison of the KT and TC phasor 
recovery techniques. Each criterion tested the reconstructed image’s resolution 
produced by both algorithms over a range of values for a specific imaging parameter. 
A baseline of imaging parameters (Table A.I) was established. Typically one 
parameter was varied within an individual criterion. Each criterion used a range of 
two to four imaging parameter values. The criteria included reconstructed image 
evaluation based on the quantity of short exposure images, the short exposure image 
photon count, and the amount of turbulence. Additionally, the effects of short 
exposure image tilt on reconstruction as well as centroiding in image space to remove 
it, were weighed. Further, the effect of offset value on cross-spectrum and 
bispectrum estimates and object size on image resolution, were rated. 

The techniques were judged in terms of image resolution, phasor spectrum 
SNR, and phase error. Each comparison included both the KT and TC image 
reconstructions. The evaluation included two-dimensional graphs of the SNR values 
for both the KT and TC image reconstruction as well as their phase error values. 
The reconstructed images appear with normalized intensities on two-dimensional 
contour plots with hachure marks indicating the direction of minima. The SNR and 
phase error values for each comparison were plotted against spatial firequent^. The 
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maximum frequency value was the last point at which the SNR had a value of one 
or greater. The reconstructed image plots were visually compared to the object’s true 
image. By contrast, their SNR and phase error values were compared to each other. 

B. RECOVERY TECHNIQUE COMPARISONS 

Appendix A maintains the results of the recoveiy technique comparisons within 
the seven criteria. Table A.I delineates the imaging parameters involved in the 
simulation process and whether these values are variable over these criteria. Figures 
A.1 and A.2 represent the true representations of the asteroid and binary star 
respectively. These figures show the results of the actual computer generated objects 
without turbulence corruption, filtering, or modification resulting from the imaging 
system including aperture effects, which all other figures include. They were the 
reference figures for the reconstructed images and the phase error calculations. 
Figures A.3 through A22 are plots of single short exposure images of the asteroid, 
the binary star and the star that show the effects of varying coherence lengths and 
photon counts. They show the level of distortion the objects realize in the imaging 
process. Figures A23 through A25 are plots of long exposure images that consisted 
of an average of 100 short exposure images with tilt. The inclusion of these images 
provides an appreciation for the necessity of image reconstruction to acquire an 
image that more closely resembles the truth. 
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1. Short Exposure Image Quantity 

Image resolution increases with the quantity of short exposure images used 
in the reconstruction process. Varying the quantity of images with no tilt present 
provided the first comparison criterion for the recovery techniques. For the fom 
comparisons, the values of Nf were 400 (baseline), 25,100, and 1600 short exposure 
images. All other baseline parameters remained unchanged. The relevant plots and 
graphs are Figures A.26 through A.41. In all cases, the visual image quality of the 
TC recovered images was superior to those recovered by the KT process. This image 
quality distinction was especially noticeable at the lower Nf values of 25 and 100. As 
the value of Nf increased, the distinction decreased to the point where it was only 
slightly noticeable at the Nf of 1600. Anafysis of the phasor spectrum SNR graphs 
showed, in general, that the TC SNR curves were offset toward approximately ten to 
20 percent greater SNR values than those of the KT SNR curves. Further, the phase 
error graphs showed the error ciirves resulting from the TC method were offset 
toward approximately five to 15 percent lesser error than those of the KT method. 
These two series of graphs confirmed that in all cases, especially at low Nf, the TC 
technique outperformed the KT technique. 

2. Photon Count 

Image resolution increases with the amount of short exposure image 
photons present. Varying the quantity of photons in the short exposure images with 
no tilt present provided the second comparison criterion for the recovery techniques. 
For the four comparisons, the photon count values were 10^ (baseline), 10^, 10^, and 
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10^ photons where all other baseline parameters remained unchanged. Figures A.26 
through A.29 and A.42 through A.S3 are applicable for this criterion. The visual 
quality of the TC and KT reconstructed images for photon counts of 10^ and 10^ was 
almost identical. For photon counts of 10^ and 10^, the KT technique produced 
slightly better reconstructed images. This photon count distinction of recovery 
technique performance was confirmed by both the SNR and phase error graphs. For 
low photon count short «posure images, the KT technique SNR curves were ofEset 
toward approximately ten to 20 percent greater SNR values than those of the TC 
technique, and the phase error curves were general^ offset toward five to 25 percent 
lesser error values. Therefore, for low photon count image recovery without tilt, the 
KT phasor recovery technique provides slight^ better resolution. 

3. Turbulence Magnitude and Offset Value 

As the amount of turbulence increases, the resolution of the reconstructed 
image decreases. Short exposure images with no tilt of coherence lengths 0.103 and 
0.206 meters (baseline), provided the third and fourth comparison criteria for the 
recovery techniques. For the turbulence magnitude criterion, the 0.103 meter 
coherence length images required an offset value of two, to ensure inclusion of all 
spatial frequencies. For consistency, the 0.206 meter coherence length images used 
the same value. Comparison between offset values of two and four satisfied the 
offset criterion for TC and KT reconstructed images. Ofifeet values effect the 
quantity of cross-spectrum or bispectrum averaging and the reconstructed image 
quality increases with increasing offset value. All other baseline parameters remained 
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unchanged. Figures A.54 through A.61 showed the comparisons of the turbulence 
criterion and Figures A.26 through A-29 and A.58 through A61 showed the 
comparisons of the offset criterion. For the coherence length case, the KT technique 
provided a slightly better image with greater turbulence. The TC SNR curve was 
offset toward approximately five to ten percent greater SNR values relative to the KT 
curve. However, the KT phase error curve was offset toward approximately five 
percent lesser error values providing the more resolved image. For the offset value 
case, the TC technique produced a slightly better image, though the SNR and phase 
error curves were almost coincident. This apparent contradiction arose from TC 
techniques having more fi-equent^ values above the unity SNR value, thereby 
providing higher frequencies for the filtering process. 

4. Tilt 

The resolution of the reconstructed image declines with the inclusion of 
tilt in the short exposure images. The previous criterion comparisons were conducted 
without the presence of tilt. The addition of tilt provides a more realistic comparison 
of the phasor recovery techniques as tilt is always present in true short exposure 
images. Varying the photon count of the short exposure images with tilt present 
provided the fifth criterion for comparison of the recovery techniques. For the three 
comparisons, the photon count values were 10^, 10*, and 10^ photons and all other 
baseline parameters remained unchanged. Figures A.62 through A.77 were 
applicable for this criterion. In all photon count cases, the TC recovered image were 
superior. The KT SNR curves were generalty offset toward 30 to 50 percent lesser 
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SNR values and dropped below the unity value at lower spatial frequencies than the 
TC curves. Hence, fewer high frequency values were included in the filtering process, 
producing a poorer resolution. The phase error curves for the TC recovered images 
were ofEset toward approximately 20 to SO percent lesser error values than for those 
recovered using the KT process except at the 10^ photon count. At this value, 
however, the KT SNR curve was offret toward much lower SNR values. 
Consequently, the TC techmque was found to be superior when tilt was included in 
the short exposure images. 

5. Centroiding 

Centroiding the short exposure images prior to image reconstruction 
enhances both the TC and KT recovery techniques for high photon counts. Short 
exposure image centroiding provided the sixth criterion for comparison of the 
recovery techniques by again varying the photon count of the short exposure images. 
With tilt present, the images were centroided prior to reconstruction. For the three 
comparisons, the photon count values were 10^, 10^, and 10^ photons and all other 
baseline parameters remained unchanged. Additionally, a comparison with 10^ 
photons and 1600 short exposures with and without tUt tested the effects of 
centroiding at high Nf values. Figures A78 through A93 were relevant for this 
criterion. Centroiding offered a two to five percent improvement to both recovery 
methods with higher photon count. At low photon counts such as 10^ photons, 
centroiding actually offset the phase error curves toward higher error values. For the 
TC method, centroiding should have no effect or the effect should disappear with 
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large enough Nf values since the method is shift invariant. However, at 10* photons 
and with 1600 short exposure images, a minor improvement in phase error occurred 
with centroiding. Centroiding did not improve the reconstructed image to thr* point 
of those recovered images having no tilt, nor did centroiding bring the KT method 
results in line with that of the TC method. However, a minor improvement in phase 
error occurred for both methods. 

6. Point Objects 

The resolution of the reconstructed image depends upon its size and 
detail. Resolution of an object such as a binary star occurs more easily because of 
the requirement for less short exposure image photons. A binary star with one star 
having twice the intensity as its counterpart, provided the seventh comparison 
criterion for the phasor recovery techniques. The binary star was corrupted by 
turbulence of coherence length 0.103 and 0.206 meters and with short exposure image 
photon counts of 10^ and 10* photons. Figures A.94 through A. 109 are germane. 
The two techniques produced identical results at the higher photon count and lower 
turbulence values. With lower photon count and greater turbulence, the TC 
technique provided greater resolution of the stars. Both the SNR and phase error 
curves supported this fact. 

7. Recovery Technique Comparison Findings 

For the eventual use of real objects in future applications of the recovery 
techniques where tilt is inherent in short exposure images, triple-correlation was the 
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superior reconstruction technique. The triple-correlation technique exceeded the 
Knox-Thompson technique by far in the comparisons where tilt was a concern. The 
shift invariance of the triple-correlation technique provided the ability to resolve low 
light-level objects where the Knox-Thompson technique failed. Though the Knox- 
Thompson technique required eight percent less time for image reconstruction, the 
resolution improvement obtained by the TC approach outweighed the computational 
time efficiency of the competing technique. 
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VI. CONCLUSION 


A. OVERVIEW 

This thesis compared the Knox-Thompson and triple-correlation phasor 
recovery techniques. Since speckle imaging involves extensive calculations, 
comparison of the two techniques required a powerful computer. The simulation 
produced an object and a phase screen for each short exposure image. The 
diffraction-limited information in these images allowed reconstruction of the object’s 
power and phasor spectra. Combining these spectra produced the reconstructed 
image. Image reconstruction comparison under seven imaging criteria permitted the 
ability to determine the superior technique. The triple-correlation technique provided 
the best overaU image resolution. This judgement stems from its superiority with 
regard to realistic short exposure images which included tilt. 

B. OPTIMUM IMAGE RECOVERY APPROACH 

This thesis found that, of the two phasor recovery techniques compared, the 
triple-correlation technique was the optimum approach for real short exposure image 
recovery. From the shift invariance of the triple-correlation technique, attainment 
of 20 to 50 percent less azimuthal phase error values occurred when compared with 
the Knox-Thompson technique. As a result, use of the triple-correlation phasor 
recovery technique is essential. Specifically, removal of the noise bias compensates 
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for the random photon noise that is intrinsic to real images. The use of phasor vice 
phase recovery is key to avoid the phase 2ir wrap-around problem. The weighting 
of the phasors determined from their cross-spectrum or bispectrum estimates by the 
least squares estimation approach, is critical. The recursion method used to provide 
the phasors is not imperative, and other techniques may be used, such as the method 
of least squares or the method of steepest decent, not discussed in this thesis. With 
regard to imaging parameters, the triple-correlation technique allowed larger 
ma ximum o&et values than the Knox-Thompson technique; D/X instead of r^/X. 
This provides maximum estimate averaging. Based on the results with centroiding, 
it is helpful for relatively high light level short exposure images. Finally, peak 
recovered image resolution requires the maximum amount of short raposure images 
practicable. 

C FURTHER STUDY 

This thesis provided simulated results for comparison of the two recoveiy 
techniques. Application of the triple-correlation technique to actual, turbulence- 
degraded images provides an avenue of research. Development and testing of other 
methods of extracting the phasor spectrum beyond the recursion method demands 
analysis. Reconstructed Fourier spectrum filtering processes beyond the simple 
rectangular low pass filtering approach used herein require exploration. 
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APPENDIX A. PLOTS AND GRAPHS 


The following table, plots, and graphs were referenced in the text. 
Figure AJ Imaging Parameters. 


Imaging Parameter Parameter Parameter Parameter 

Abbrev. Baseline Value 

Short Exposure Nf 400 Variable 

Image Quantity 

Coherence Length rg 0.206 (m) Variable 

Short Exposure Np 10* Variable 

Image Photon 

Quantity 

Offset Value OS 4 (pixels) Variable 

Telescope Primary D 1.6 (m) Constant 

Diameter 

Telescope Secondary D, 0.33 (m) Constant 

Diameter 

Aperture Radius a 31 (pixels) Constant 

Light Wavelength k 5,5 x 10'^ Constant 

(m) 

Cutoff Frequency f^ 68.3 (1/arcsec) Constant 

Random Number s 123456789 Constant 

Seed 
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Figure A3 Asteroid Short Exposure Image: 

ro = 0.206, Np = 10^. 



arc—seconds 

Figure A.4 Asteroid Short Exposure Image, 

ro = 0.206, Np = 10*. 
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Figure AS Asteroid Short Exposure Image, 


ro = 0.206, Np = 10*. 



arc—seconds 

Figure A.6 Asteroid Short Exposure Image, 

ro = 0.206, Np = 10*. 


54 




0.55 


m 

rt 

G 

o 

O _ 


- 0.55 


- 1.10 



- 1.10 - 0.55 - 0.00 0.55 

arc-seconds 

Figure A.7 Asteroid Short Exposure Image, 
To = 0.103, Np = 10*. 
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Figiire Asteroid Short Exposure Image, 
ro = 0.103, N_ * 10*. 
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Figure A.10 Asteroid Short Exposure Image, 
ro = 0.103, Np = lO*. 
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Figure A.11 Binary Star Short Exposure Image, 
ro = 0.206, Np = lO*. 



Figure A.12 Binary Star Short Exposure Image, 

ro = 0.206, Np = 10^, 


57 
























J_I_I_I_I 

-1.10 -0.55 -0.00 0.55 

arc-seconds 


Figure A.18 Star Short Exposure Image, 
ro = 0.206, Np = lO^. 
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Figure A.19 Star Short Exposure Image, 
ro = 0.103, Np = 10*. 




Figure A.20 Star Short Exposure Image, 
ro = 0.103, Np = 10^. 
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Figure K22 Star Short Exposure Image, 
ro = 0.103, Np = lO^. 
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Figure A^<> Baseline KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 400, OS = 4. 
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Figure A.27 Baseline TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 400, OS = 4. 
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Figure A^8 Asteroid Baseline Phasor Spectrum SNR. 



Spatial Frequency (l/arcsec) 

Figure AJ29 Asteroid Baseline Azimuthal RMS Phase Error. 
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Figure A30 KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 25, OS = 4. 



Figure A31 TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 25, OS = 4. 
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Figure A32 Asteroid Phaser Spectrum SNR, Nf = 25. 
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Figure A33 Asteroid Azimuthal RMS Phase Error, Nf = 25. 
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Figure AM KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 100, OS = 4. 
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Figure AJ5 TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, Nj = 100, OS = 4. 
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Figure AJ6 Asteroid Phasor Spectrum SNR, Nf = 100. 



Figure A37 Asteroid Azimuthal RMS Phase Error, Nf = 100. 
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Figure A38 KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 1600, OS = 4. 
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Figure AJ9 TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10®, N, = 1600, OS = 4. 
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Figure A.40 Asteroid Phasor Spectrum SNR, Nj = 1600. 



Figure A.41 Asteroid Azimuthal RMS Phase Error, Nf = 1600. 


72 







o.ss 

m 

•O 

C 

O 

V -0.00 
CO 

I 

o 

u, 

«a 

- 0.66 


- 1.10 

- 1.10 - 0.66 - 0.00 0.66 
arc-seconds 

Figure A.42 KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10®, N, = 400, OS = 4. 
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gun A.43 TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10^, N, = 400, OS = 4. 
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Figure A.44 Asteroid Phasor Spectrum SNR, Np = 10*^. 



Figure A.45 Asteroid Azimuthal RMS Phase Error, Np = 10^. 
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Figure A.46 KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10^, N, = 400, OS = 4. 
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Figure A.47 TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10^, N, = 400, OS = 4. 
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Figure A.48 Asteroid Phaser Spectrum SNR, Np = 10^. 



Figure A.49 Asteroid Azimuthal RMS Phase Error, Np = 10^. 
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Figure AJO KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = lO*, N, = 400, OS = 4. 



Figure A.51 TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10^, N^ = 400, OS = 4. 


77 







RMS Phase Error (radians) Phaser Spectrum SNR 



Figure AJ2 Asteroid Phaser Spectrum SNR, Np = 1(P. 



Figure A.53 Asteroid Azimuthal RMS Phase Error, Np = 10^. 
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Figure A.54 KT Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 400, OS = 2. 
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Figure A.55 TC Recovered Asteroid, No Tilt, 
ro = 0.206, Np = 10*, N, = 400, OS = 2. 
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Figure A.56 Asteroid Phaser Spectrum SNR, OS = 2. 



Figure A.57 Asteroid Azimuthal RMS Phase Error, OS = 2. 
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Figure A^8 KT Recovered Asteroid, No Tilt, 


ro = 0.103, Np = 10*, Nf = 400, OS = 2. 
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Figure A.59 TC Recovered Asteroid, No Tilt, 
ro = 0.103, Np = 10*, Nf = 400, OS = 2. 
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Figure A.60 Asteroid Phaser Spectrum SNR, 
ro = 0.013, OS = 2. 



Spatial Frequency (1/arcsec) 

Figure A.61 Asteroid Azimuthal RMS Phase Error, 
ro = 0.103, OS = 2. 
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Figure A.62 KT Recovered Asteroid, 
With Tilt, To = 0.206, N = 10*, 

N, = 400, OS = 4. 
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Figure A.63 TC Recovered Asteroid, 
With Tilt, ro = 0.206, N = 10*, 

Ni = 400, OS = 4. 
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Figure A.64 Asteroid Phaser Spectrum SNR, Tilt. 
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Figure A.65 Asteroid Azimuthal RMS Phase Error, Tilt 
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Figure A.68 Asteroid Phaser Spectrum SNR, 
Tilt, Np 10^. 
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Figure A.69 Asteroid Azimuthal RMS Phase Error, 
Tilt, Np * 10^. 
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Figure A.70 KT Recovered Asteroid, 
With Tilt, ro = 0.206, N = lO*, 

Nf = 400, OS = 4. 



Figure A.71 TC Recovered Asteroid, 
With Tilt, ro = 0.206, Np = 1(>*, 

Nf = 400, OS = 4. 
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Figure A.72 Asteroid Phaser Spectrum SNR, 
Tilt, Np = icy*. 



Figure A.73 Asteroid Azimuthal RMS Phase Error, 
Tilt, Np = 10^. 
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Figure A.74 KT Recovered Asteroid, 
With Tilt, ro = 0.206, N = 10*, 

Nf = 1600, OS = 4. 
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Figure A.75 TC Recovered Asteroid, 
With Tilt, ro = 0.206, N = 10*, 

Nt = 1600, OS = 4. 
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Figure A.78 KT Recovered Asteroid, 
With Tilt, Centroided, ro = 0.206, 
Np = 10*, N, = 1600, OS = 4. 
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Figure A.79 TC Recovered Asteroid, 
With Tilt, Centroided, r© = 0.206, 
Np = 10*, N, = 1600, OS = 4. 
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Figure A.80 Asteroid Phaser Spectrum SNR, 
Tilt, Centroiding, Nj = 1600. 



Spatial Frequency (l/arcsec) 

Figure A.81 Asteroid Azimuthal RMS Phase Error, Tilt, 
Centroiding, Nf = 1600. 
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Figure A.82 KT Recovered Asteroid, 
With Tilt, Centroided, r® = 0.206, 
Np = 10*, Nf = 400, OS = 4. 




-i.io I-1-1-1_I_I_I_t 

- 1.10 - 0.66 - 0.00 0.66 
arc—seconds 

Figure A.83 TC Recovered Asteroid, 
With Tilt, Centroided, ro = 0.206, 

Np = 10*, Nf = 400, OS = 4. 
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Figure A.84 Asteroid Phaser Spectrum SNR, Tilt, Centroiding. 



Figure A.85 Asteroid Azimuthal RMS Phase Error, Tilt, Centroiding. 
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Figure A.86 KT Recovered Asteroid, With Tilt, 
Centroided, r® = 0.206, N. = 10^, 

Nf = 400, OS = 4. 



arc-seconds 

Figure A.87 TC Recovered Asteroid, With Tilt, 
Centroided, r® = 0.206, N, = 10*, 

Nf = 400, OS = 4. 
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Phasor Spectrum SNR 










Figure A-90 KT Recovered Asteroid, With Tilt, 
Centroided, r^, = 0.M6, N. = lO*, 

N, = 400, OS = 4. 



Figure A.91 TC Recovered Asteroid, With Tilt, 
Centroided, ro = 0.206, N- = 10^, 

Nj = 400, OS = 4. 
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RMS Phase Error (radians) Phaser Spectrum SNR 



Figure A.96 Binaiy Star Phaser Spectrum SNR, Np = 1(P. 



Figure A.97 Binaiy Star Azimuthal RMS Phase Error, Np = 10^. 
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1.55 



Figure A.98 KT Recovered Binary Star, No Tilt, 
ro = 0.206, Np = lO^, N, = 400, OS = 4. 



Figure A.99 TC Recovered Binary Star, No Tilt, 
ro = 0.206, Np = lO^, N, = 400, OS = 4. 
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RMS Phase Error (radians) Phasor Spectrum SNR 



Figure A.100 Binaiy Star Phasor Spectrum SNR, Np = 10^. 



Figure A.101 Binary Star Azimuthal RMS Phase Error, Np = 10^. 
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Figure A.102 KT Recovered Binaiy Star, No Tilt, 
To = 0.103, Np = IO3, Nf = 400, OS = 4. 



Figure A.103 TC Recovered Binaiy Star, No Tilt, 
ro = 0.103, Np = IO3, N, = 400, OS = 4. 
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Figure A.104 Binary Star Phaser Spectrum SNR, r^ = 0.103, Np = 10 



Figure A.105 Binary Star Azimuthal RMS Phase Error, 







arc-seconds 


Figure A.106 KT Recovered Binary Star, No Tilt, 
ro = 0.103, Np = 10^, Nf = 400, OS = 4. 



arc-seconds 

Figure A.107 TC Recovered Binary Star, No Tflt, 
ro = 0.103, Np = 10^, Nf = 400, OS = 4. 
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Figure A.108 Binary Star Phaser Spectrum SNR, r^ = 0.103, Np = lO' 



Figure A.109 Binary Star Azimuthal RMS Phase Error, 







APPENDK B. KNOX-THOMPSON MAIN PROGRAM 


c THIS PROGRAM CREATES ONE OF THREE IMAGES: A STAR, A 
c BINARY STAR AND AN ASTEROID. IT THEN DEGRADES THE 
c IMAGE BY SIMULATING ATMOSPHERIC CONDITIONS AND PHOTON 

c NOISE, AND THEN RECONSTRUCTS THE IMAGE USING THE 
c KNOX-THOMPSON ALGORITHM. 


c 

c 

c 

c 

c 

c 

c 

c 


AUTHOR: LT JAMES M. LACKEMACHER 

COMPL DATE: 26 OCTOBER 1990 

REASON: COMPLETE REQUIREMENTS FOR A MASTERS 

DEGREE IN PHYSICS. 

GOAL: SIMULATE OBJECT, DEGRADE OBJECT, 

RECONSTRUCT OBJECT USING KNOX- 
THOMPSON AND TRIPLE-CORRELATION 
METHODS, FILTER AND COMPARE. 


PROGRAM KNOXTHOMPSON 


c MAIN PROGRAM COMPLEX VARIABLE LIST 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


F n DIM ARRAY USED IN THE FOURIER TRANSFORM 

GAUSSIAN n X n DIM ARRAY THAT REPRESENTS THE GAUSSIAN 
PORTION OF THE ASTEROID 
I n X n DIM ARRAY THAT REPRESENTS THE 

DEGRADED IMAGE 

nCT nxnx5x9 ARRAY THAT REPRESENTS THE CROSS¬ 

SPECTRUM OF THE IMAGE 

IS n X n DIM ARRAY THAT REPRESENTS THE 

DEGRADED POINT SOURCE 

ISKT n X n DIM ARRAY THAT REPRESENTS THE MODULUS 

SQUARED OF THE POINT SOURCE 


107 





c 

c 

c 

c 

c 

c 

c 


c 


OKT 

OBJDATA 

KTPHASOR 

PUPIL 

TEMPDATA 


n X n DIM ARRAY THAT REPRESENTS THE OBJECT 
SPECTRUM OF THE RECONSTRUCTED IMAGE 
n X n DIM ARRAY THAT REPRESENTS THE OBJECT 
n X n DIM ARRAY THAT REPRESENTS THE PHASOR OF 
THE OBJECT IN THE RECONSTRUCHON PROCESS 
n X n DIM ARRAY THAT REPRESENTS THE PUPIL 
PORTION OF THE ASTEROID 
n X n DIM ARRAY THAT IS USED AS A TEMPORARY 
ARRAY 


c 


MAIN PROGRAM REAL VARIABLE UST 


c KTsnr 
c 

c mod 
c 

c rsnr 
c 

c xvarrkt 

c 

c 

c xvarikt 

c 

c 


n X n DIM ARRAY THAT REPRESENTS THE SNR OF 
EACH PHASOR 

n X n DIM ARRAY THAT REPRESENTS THE MODULUS 
OF THE RECONSTRUCTED IMAGE 
n/2 DIM ARRAY THAT REPRESENTS THE SNR AS A 
FUNCTION OF RADIUS 

n X n X 5 X 9 DIM ARRAY THAT REPRESENTS THE 
REAL PART OF THE VARIANCE OF THE CROSS- 
SPECTRUM 

n X n X 5 X 9 DIM ARRAY THAT REPRESENTS THE 
IMAGINARY PART OF THE VARIANCE OF THE CROSS- 
SPECTRUM 


c MAIN PROGRAM INTEGER VARIABLE UST 


c fwd 
c icounter 
c 

c inv 
c In 
c o^et 
c 
c 

c mseed 
c 


VALUE OF 1 FOR FORWARD FFT 
COUNTER THAT COUNTS THE NUMBER OF 
SNAPSHOTS 

VALUE OF -1 FOR INVERSE FFT 
2^^^ FOR USE WITH FFT SUBROUTINE 
VARIABLE THAT REPRESENTS THE NUMBER OF 
PIXELS THAT ARE AVERAGED IN THE KNOX- 
THOMPSON PROCESS 

VARIABLE USED TO ENSURE ONLY ONE PASS OF 
INITIAL PART OF PHSUB SUBROUTINE 
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c n 
c nframes 
c nphoton 
c nyquist 
c 
c 
c 


DIMENSION OF ONE SIDE OF 2-DIM ARRAY 
TOTAL NUMBER OF SHORT EXPOSURE SNAPSHOTS 
TOTAL NUMBER OF PHOTONS IN SNAPSHOT 
EQUAL TO THE TELESCOPE PUPIL FUNCTION RADIUS 
DERIVED FROM THE RELATION: 
nyquist = (telescope diameter x 
number of pixels per rO)/rO 


MAIN PROGRAM 


PARAMETER(ii=64,ln=6,fwd- l,inv=-l) 

COMPLEX* 16 OBJDATA(n,n), TEMPDATA(n,n), F(n), 

+ PUPIL(n,n), STARDATA(n,n). GAUSSIAN(n,n), 

+ KTPHASOR(n,n), IKT(n,n,5,9), ISKT(n,n), 

+ I(n,n), IS(n,n), OKT(n,n) 

REAL*8 mod(n,n), xvarrkt(n,n,5,9), xvarilct(n,n,5,9), 

+ KTsnr(n,n), rsnr(n/2) 

INTEGER offset 

CHARACTER* 16 filel, file2 
CHARACTER*! cent 

c INITIALIZE MSEED TO ALLOW ONLY ONE PASS THROUGH FIRST 
c PART OF PHSUB 

mseed = 1 

c INITIALIZE PROGRAM READING REQUIRED VARIABLES AND 
c CREATE THE DESIRED OBJECT 

CALLInitialize(OBJDATA,GAUSSIAN,PUPIL,F,TEMPDATA, 

+ filel,file2,fwd,inv,In,nframes,nphoton,nyquist, 

+ offset,cent,n) 

c TRANSFORM THE OBJECT TO FREQUENCY SPACE 
CALL FFT2D(OBJDATA,TEMPDATA,F,ln,fwd,n) 
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c NOR^MLIZE ALL OF THE ARRAY ELEMENTS TO THE VALUE OF 
c THE NUMBER OF PHOTONS WHERE THE DC TERM EQUALS THE 
c NUMBER OF PHOTONS OF THE IMAGE 

CALL Photons(OBJDATA,nphoton^) 

c CREATE THE POINT SOURCE 

CALL Star(STARDATA,n) 

c TRANSFORM THE POINT SOURCE TO FREQUENCY SPACE 

CALL FFT2D(STARDATA,TEMPDATA,F4n,fwd,n) 

c NORMALIZE ALL OF THE ARRAY ELEMENTS TO THE VALUE OF 
c THE NUMBER OF PHOTONS WHERE THE DC TERM EQUALS THE 
c NUMBER OF PHOTONS OF THE IMAGE 

CALL Photons(STARDATA,nphoton,n) 

c COMMENCE THE LOOP THAT COUNTS THE SNAPSHOTS 

DO 10 icounter « 1, nframes 

c DEGRADE THE IMAGE WITH THE TELESCOPE, THE ATMOSPHERE, 
c AND THE PHOTON NOISE 

CALL PHSUB(I,F,TEMPDATA,OBJDATA,icounter,nframes, 

+ nyquist,fwd,inv,ln,n,inseed) 

c CENTROID THE DEGRADED IMAGE ONLY SINCE ONLY PHASE IS 
c EFFECTED BY CENTROIDING IF CENTROIDING IS DESIRED. 

IF ((ccnt.EQ.’Y’).OR.(cent.EQ.y)) THEN 
CALL Centroid(I,TEMPDATA,n) 

ENDIF 

c TRANSFORM THE DEGRADED IMAGE TO FREQUENCY SPACE 
CALL FFT2D(I,TEMPDATA,F,ln,fwd,n) 


no 






c DEGRADE THE POINT SOURCE WITH THE TELESCOPE, THE 
c ATMOSPHERE, AND THE PHOTON NOISE 

CALLPHSUB(IS,F,TEMPDATA,STARDATA,icounter,nframes, 

+ nyquist,fwd,inv,ln,n,inseed) 

c TRANSFORM THE DEGRADED IMAGE TO FREQUENCY SPACE 

CALL FFT2D(IS,TEMPDATA,F,In,fwd,n) 

c DETERMINE THE MODULUS OF THE DEGRADED IMAGE AND 
c POINT SOURCE 

CALL Modulus(I,IS,IKT,ISKT,mocl,icounter, 

+ nfirames,nphoton,n) 

c RECONSTRUCT THE OBJECT FROM THE DEGRADED IMAGE AND 
c POINT SOURCE 

CALL KTrecon(KTPHASOR,I,IKT,KTsnr,offset,xvarrkt, 

+ xvarikt,icounter,nfranies,nyquist,n) 

WRITE(*,*)icounter,TRAMES COMPLETED’ 

c END THE LOOP 

10 CONTINUE 

c CALCULATE THE AVERAGE SNR AS A FUNCTION OF RADIUS 

CALL SNRcalc(KTsnr,rsnr,nyquist,n) 

c COMBINE THE MODULUS WITH THE PHASOR 

CALL Combine(OKT,raod,KTPHASOR,n) 

c WRITE RECONSTRUCTED PHASE AND POWER SPECTRUM TO A 
c FILE 

CALL Writefile(OKT,rsnr,file l,file2,nyquist,n) 

STOP 

END 
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APPENDIX C TRIPLE-CORRELATION MAIN PROGRAM 


c THIS PROGRAM CREATES ONE OF THREE IMAGES: A STAR, A 
c BINARY STAR AND AN ASTEROID. IT THEN DEGRADES THE 

c IMAGE BY SIMULATING ATMOSPHERIC CONDITIONS AND PHOTON 

c NOISE, AND THEN RECONSTRUCTS THE IMAGE USING THE TRIPLE- 

c CORRELATION ALGORITHM. 


c 

c 

c 

c 

c 

c 

c 

c 


AUTHOR: LT JAMES M. LACKEMACHER 

COMPL DATE: 26 OCTOBER 1990 

REASON: COMPLETE REQUIREMENTS FOR A MASTERS 

DEGREE IN PHYSICS. 

GOAL: SIMULATE OBJECT, DEGRADE OBJECT, 

RECONSTRUCT OBJECT USING KNOX- 
THOMPSON AND TRIPLE-CORRELATION 
METHODS, FILTER AND COMPARE. 


PROGRAM TRIPLECORR 


c MAIN PROGRAM COMPLEX VARIABLE LIST 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BSPHASOR 

F 

GAUSSIAN 

I 

IBS 

IDBS 


n X n DIM ARRAY THAT REPRESENTS THE PHASOR OF 
THE OBJECT IN THE RECONSTRUCTION PROCESS 
n DIM ARRAY USED IN THE FOURIER TRANSFORM 
n X n DIM ARRAY THAT REPRESENTS THE GAUSSIAN 
PORTION OF THE ASTEROID 

n X n DIM ARRAY THAT REPRESENTS THE DEGRADED 
IMAGE 

n X n X 5 X 9 ARRAY THAT REPRESENTS THE 
BISPECTRUM OF THE IMAGE 
5x9 ARRAY THAT REPRESENTS THE DEVIATION 
FROM THE DC TERM 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IS 

ISBS 

OBS 

OBJDATA 

PUPIL 

TEMPDATA 


n X n DIM ARRAY THAT REPRESENTS THE DEGRADED 
POINT SOURCE 

n X n DIM ARRAY THAT REPRESENTS THE MODULUS 
SQUARED OF THE POINT SOURCE 
n X n DIM ARRAY THAT REPRESENTS THE OBJECT 
SPECTRUM OF THE RECONSTRUCTED IMAGE 
n X n DIM ARRAY THAT REPRESENTS THE OBJECT 
n X n DIM ARRAY THAT REPRESENTS THE PUPIL 
PORTION OF THE ASTEROID 
n X n DIM ARRAY THAT IS USED AS A TEMPORARY 
ARRAY IN THE FOURIER TRANSFORM 


c 


MAIN PROGRAM REAL VARIABLE LIST 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BSsnr 

mod 

rsnr 

xvarrbs 

xvaribs 


n X n DIM ARRAY THAT REPRESENTS THE SNR OF 
EACH PHASOR 

n X n DIM ARRAY THAT REPRESENTS THE 
MODULUS OF THE RECONSTRUCTED IMAGE 
n/2 DIM ARRAY THAT REPRESENTS THE SNR AS A 
FUNCTION OF RADIUS 

n X n X 5 X 9 DIM ARRAY THAT REPRESENTS THE 
REAL PART OF THE VARIANCE OF THE BISPECTRUM 
n X n X 5 X 9 DIM ARRAY THAT REPRESENTS THE 
IMAGINARY PART OF THE VARIANCE OF THE 
BISPECTRUM 


c MAIN PROGRAM INTEGER VARIABLE LIST 


c fwd 

c icounter 

c 

c inv 

c In 

c offset 


VALUE OF 1 FOR FORWARD FFT 
COUNTER THAT COUNTS THE NUMBER OF 
SNAPSHOTS 

VALUE OF -1 FOR INVERSE FFT 
2"'In FOR USE WITH FFT SUBROUTINE 
VARIABLE THAT REPRESENTS THE NUMBER OF 
PIXELS THAT ARE AVERAGED IN TFE TRIPLE¬ 
CORRELATION PROCESS 
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+ + 


c mseed 
c 

c n 
c nframes 
c nphoton 
c nyquist 
c 


VARIABLE USED TO ENSURE ONLY ONE PASS OF 
INITIAL PART OF PHSUB SUBROUTINE 
DIMENSION OF ONE SIDE OF 2-DIM ARRAY 
TOTAL NUMBER OF SHORT EXPOSURE SNAPSHOTS 
TOTAL NUMBER OF PHOTONS IN SNAPSHOT 
EQUAL TO THE TELESCOPEPUPIL FUNCTION RADIUS 
DERIVED FROM THE RELATION: 


c nyquist = (telescope diameter x 

c number of pixels per rO)/iO 


MAIN PROGRAM 


PARAMETER(n=64,ln=6,fwd=l,inv=-l) 


COMPLEX*16 OBJDATA(n,n), TEMPDATA(n,n), F(n), 

+ PUPIL(n,n), I(n,n), STARDATA(n,n), GAUSSIAN(n,n), 

+ BSPHASOR(n,n), IBS(n4i,5,9), ISBS(n,n), IS(n,n), 

+ OBS(n,n), IDBS(5,9) 

REAL*8 mod(n,n), xvarrbs(n^,5,9), xvaribs(n,n,5,9), 

+ BSsnr(n,n), rsnr(n/2) 

INTEGER offset 

CHARACTER*16 fflel, file2 
CHARACTER* 1 cent 

c DSrmALIZE MSEED TO ALLOW ONLY ONE PASS THROUGH FIRST 
c PART OF PHSUB 

mseed = 1 

c INITIALIZE PROGRAM READING REQUIRED VARIABLES AND 
c CREATE THE DESIRED OBJECT 

CALLInitialize(OBJDAT>^GAUSSIAN,PUPIL,F,TEMPDATA, 
filel,file2,fwd,inv,ln,inframes, 
nphoton,nyquist,offseLcent,n) 
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c TRANSFORM THE OBJECT TO FREQUENCY SPACE 

CALL FFT2D(OBJDATA,TEMPDATA,F,ln^d,n) 

c NORMALIZE ALL OF THE ARRAY ELEMENTS TO THE VALUE OF 
c THE NUMBER OF PHOTONS WHERE THE DC TERM EQUALS THE 
c NUMBER OF PHOTONS OF THE IMAGE 

CALL Photons(OBJDATA,nphotoiMi) 

c CREATE THE POINT SOURCE 

CALL Star(STARDATA,n) 

c TRANSFORM THE POINT SOURCE TO FREQUENCY SPACE 

CALL FFT2D(STARDATA,TEMPDATA,F4n,fwd,n) 

c NORMALIZE ALL OF THE ARRAY ELEMENTS TO THE VALUE OF 
c THE NUMBER OF PHOTONS WHERE THE DC TERM EQUALS THE 
c NUMBER OF PHOTONS OF THE IMAGE 

CALL Photoiis(STARDATA,nphoton,n) 

c COMMENCE THE LOOP THAT COUNTS THE SNAPSHOTS 

DO 10 icounter = 1, nframes 

c DEGRADE THE IMAGE WITH THE TELESCOPE, THE ATMOSPHERE, 

c AND THE PHOTON NOISE 

CALL PHSUB(I,F,TEMPDATA,OBJDATA,icounter,nfirames, 

+ nyquist,fwd,inv,ln,n,mseed) 

c CENTROID THE DEGRADED IMAGE ONLY SINCE ONLY PHASE IS 
c EFFECTED BY CENTROIDING IF CENTROIDING IS DESIRED. 

IF ((cent.EQ.’Y’).OR.(cent.EQ.y)) THEN 
CALL Centroid(I,TEMPDATA,n) 

ENDIF 
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+ + 


c TRANSFORM THE DEGRADED IMAGE TO FREQUENCY SPACE 

CALL FFT2D(I,TEMPDATA,F.ln4wd,n) 

c DEGRADE THE POINT SOURCE WITH THE TELESCOPE, THE 
c ATMOSPHERE, AND THE PHOTON NOISE 

CALL PHSUB(IS,F,TEMPDATA,STARDATA,icounter,nframes, 

+ nyquist,fwd,inv,ln,n,mseed) 

c TRANSFORM THE DEGRADED IMAGE TO FREQUENCY SPACE 

CALL FFT2D(IS,TEMPDATA,F4n4wd,n) 

c DETERMINE THE MODULUS OF THE DEGRADED IMAGE AND 

c POINT SOURCE 

CALL Modulus(I,IS,IBS,ISBS 4 nod,icounter, 

+ nframes,nphoton,n) 

c RECONSTRUCT THE OBJECT FROM THE DEGRADED IMAGE AND 

c POINT SOURCE 

CALL BSrecon(BSPHASOR,I,IBS,IDBS,BSsnr,offset, 
xvarrbs,xvaribs,icounter,n£rames, 
nyquist,nphoton,n) 

WRITE(*,*)icounter,TRAMES COMPLETED’ 
c END THE LOOP 
10 CONTINUE 

c CALCULATE THE AVERAGE SNR AS A FUNCTION OF RADIUS 

CALL SNRca]c(BSsnr,rsnr,nyquist,n) 
c COMBINE THE MODULUS WITH THE PHASOR 
CALL Combine(OBS,mod,BSPHASOR,n) 
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c WRITE RECONSTRUCTED PHASE AND POWER SPECTRUM TO A 

c FILE 


CALL Writefile(OBS^nr,filel^e2,nyquistyn) 

STOP 

END 
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APPENDK D. UNIVERSAL IMAGE RECOVERY SUBROUTINES 


c THE FOLLOWING SUBROUTINES WERE GENERATED BY THE 

c THESIS AUTHOR AND ARE REQUIRED BY BOTH THE KNOX- 
c THOMPSON AND THE TRIPLE-CORRELATION PROGRAMS, 
c ADDITIONALLY, SOME SUBROUTINES ARE REQUIRED FOR 
c THE FILTRMS PROGRAM IN APPENDIX H. 


SUBROUTINE LIST 


SUBROUTINE Ast(ASTEROID,GAUSSIAN,PUPIL,F,TEMPDATA, 
+ fwd,inv4n,n) 

COMPLEX*16 GAUSSIAN(n4i), ASTEROID(n,n), PUPIL(n,n), 

+ TEMPDATA(n,n), F(ii), El, E2, E3, E4, E5, E6, E7 
maxval » 0.0 
aval = 4.0 
n2pl = 11 / 2+1 
rad = 16.0 
DO 10 i = 1, n 
DO 10 j = l,n 
X = float(j - (n2pl)) 
y = £loat((n2pl) - i) 
radius = sqrt(x**2.0 + y**2.0) 

IF (radius.LE.rad) THEN 
PUPIL(y) = (1.0,0.0) 

ELSE 

PUPIL(id) = (0.0,0.0) 

ENDIF 

IF ((abs(x).LE.aval).AND.(abs(y).LE.aval)) THEN 
GAUSSIAN(id) = DCMPLX(exp(-(x**2.0 + 

+ y**2.0)/4)) 

ELSE 

GAUSSIAN(iJ) = (0.0,0.0) 

ENDIF 

10 CONTINUE 

CALL FFT2D(GAUSSIAN,TEMPDATA,F,ln,fwd,n) 

CALL FFT2D(PUPIL,TEMPDATA,F,In,fwd,n) 
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DO 20 i = 1, n 
DO20j = l,n 

ASTEROID(ij) = GAUSSIAN(io) * PUPIL(ij) 

20 CX)NTINIJE 

CALL FFT2D(ASTEROID,TEMPDATAF,ln,mv,n) 
a = 10 
al = 2.5 
a2 = 2 
a3 = 2.5 
a4 = 2 
a5 = 3.5 
a6 = 3 
a? = 2.5 
DO 30 i = 1, n 
DO 30 j = 1, n 
X = float(j - (n2pl)) 
y = float((n2pl)- i) 
xl = x+6 
yl = y-8 
xal = xl/al 
yal = yl/al 

IF ((abs(xal).LE.aval)AND.(abs(yal).LE.aval)) 

+ THEN 

El = DCMPLX(a * exp(.(xal**2.0 + yal**2.0))) 
ELSE 

El = (0.0,0.0) 

ENDIF 
x2 = x-4 
y2 = y-6 
xa2 = x2/a2 
ya2 = y2/a2 

IF ((abs(xa2).LE.aval)AND.(abs(ya2).LE.aval)) 

+ THEN 

E2 = DCMPLX(a * exp(-(xa2**2.0 + ya2**2.0))) 
ELSE 

E2 = (0.0,0.0) 

ENDIF 
x3 = x-10 
y3 = y-4 
xa3 = x3/a3 
ya3 = y3/a3 
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IF ((abs(xa3).LE.aval)AND.(abs(ya3).LE.aval)) 

+ THEN 

E3 = DCMPLX(a * exp(-(xa3**2.0 + ya3**2.0))) 
ELSE 

E3 = (0.0,0.0) 

ENDIF 
x4 = x+6 
y4 = y+6 
xa4 « x4/a4 
ya4 « y4/a4 

IF ((abs(xa4).LE.avaI).AND.(abs(ya4)X£.aval)) 

+ THEN 

E4 = DCMPLX(a * exp(-(xa4**2.0 + ya4**2.0))) 
ELSE 

E4 = (0.0,0.0) 

ENDIF 
x5 = x+0 
y5 = y.2 
xa5 = x5/a5 
ya5 = y5/a5 

IF ((abs(xa5).LE.aval).AND.(abs(ya5).LE.aval)) 

+ THEN 

E5 = DCMPLX(a * exp(-(xa5**2.0 + ya5**2.0))) 
ELSE 

E5 = (0.0,0.0) 

ENDIF 
x6 = x+2 
y6 = y+8 
xa6 = x6/a6 
ya6 = y6/a6 

IF ((abs(xa6).L£.aval).AND.(abs(ya6).LE.ava])) 

+ THEN 

E6 = DCMPLX(a • exp(-(xa6**2.0 + ya6**2.0))) 
ELSE 

E6 = (0.0,0.0) 

ENDIF 
x7 = x-6 
y7 = y+6 
xa7 = x7/a7 
ya7 = y7/a7 
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IF ((abs(xa7).LE.aval)j\ND.(abs(ya7).LE.aval)) 

+ THEN 

E7 = DCMPLX(a • exp(-(xa7**2.0 + ya7**2.0))) 
ELSE 

E7 = (0.0,0.0) 

ENDIF 

ASTEROID(ij) = DCMPLX(DREAL(ASTEROID(ij) - 
+ ^1 + E2 + E3 + E4 + E5 + E6 + E7))) 

IF (DREAL(ASTEROID(ij)).LT.O.O) ASTEROID(ij) = 

+ (0.0,0.0) 

30 CONTINUE 
RETURN 
END 


SUBROUTINE Bistar(DATA,n) 

c THIS S/R CREATES A SIMULATED BINARY STAR WITH ONE STAR 
c LARGER THAN THE OTHER 

COMPLEX*16 DATA(n,n) 
n2pl = n/2 + 1 
DO 10 i = 1, n 
DO10j = l,n 
X * float(j - n2pl) 
y = float(n2pl - i) 

IF ((x.EQ.12.0).AND.(y.EQ.12.0)) THEN 
DATA(ij) = (2.0,0.0) 

ELSEIF ((x.EQ.-12.0)j\ND.(y.EQ.-12.0)) THEN 
DATA(ij) = (1.0,0.0) 

ELSE 

DATA(iJ) = (0.0,0.0) 

ENDIF 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE Centroid(DATA,TEMPDATA,n) 

c THIS S/R DETERMINES THE CENTROID OF THE DEGRADED IMAGE 

c AND THEN CENTROIDS THE IMAGE 

COMPLEX*16 DATA(n,n), TEMPDATA(n,n) 

REAL*8 yi, xnum, ynum, xden, yden 
n2pl = 11 / 2+1 
DO 10 i = 1, n 
DO 10 j = 1, n 
ig = dfloatQ - n2pl) 
yi = d£Ioat(ii2pl • i) 
xnum = xnum + jg * ABS(DATA(ij)) 
ynum = ynum + yi * ABS(DATA(ij)) 
xden = xden + ABS(DATA(y)) 
yden = yden + ABS(DATA(ij)) 

10 CONTINUE 

jxbar = idnint(xnum^den) 
iybar = idnint(ynum^den) 

DO 20 i = 1, n 
DO20j = l,n 
ii = i - iybar 
jj = j + jxbar 

IF((ii.GT.n).OR.(ii.LT.l).OR.(ij.GT.n).OR. 

+ (ij.LT.l)) THEN 

TEMPDATA(ij) = (0.0,0.0) 

ELSE 

TEMPDATA(ij) = DATA(iiaj) 

ENDIF 

20 CONTINUE 
DO 30 i = 1, n 
DO30j = l,n 

DATA(io) = TEMPDATA(ij) 

30 CONTINUE 
RETURN 
END 
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SUBROUTINE Combme(0,mod,PHASOR,n) 

c THIS S/R CX)MBINES THE MODULUS AND PHASOR OF THE 
c RECONSTRUCTED OBJECT 

COMPLEX*16 PHASOR(-ny2:ny2-l,-iiy2:ii/2-l), 

+ 0(-n/2:n/2-l,-n/2:n/2-l) 

REAL*8 mod(-n/2:n/2-l,-n/2:n/2-l) 
n2 = ii/2 
n2ml = ii2 -1 
DO 10 i = -n2, n2ml 
DO 10 j = -Til, n2ml 
0(y) = mod(ij) ♦ PHASOR(ij) 

10 CONTINUE 
RETURN 
END 

SUBROUTINE Complexconj(DATA,n) 

c THIS S/R IS CALLED BY FFT2D AND DETERMINES THE COMPLEX 
c CONJUGATE OF THE 2-D ARRAY IN ORDER FOR THE ARRAY TO 
c BE INVERSE FFTed. 

COMPLEX*16 DATA(n,n) 

DO 10 i = 1, n 
DO 10 j = 1, n 

DATA(id) = DCONJG(DATA(id)) 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE FFT(F,ln,n) 

c THIS S/R IS USED BY THE 2-D FFT TWICE, FIRST FFTing THE 
c ROWS THEN THE COLUMNS OF THE 2-D ARRAY. THIS S/R WAS 
c AQUIRED FROM 'DIGITAL IMAGE PROCESSING" BY GONZALEZ 
c AND WINTZ. 

COMPLEX*16 F(n),U,W,T 
REAL*8 pi, one 
one = l.OD+00 
pi = DACOS(-one) 
nv2 = n/2 
nml = n -1 

j = 1 

DO 3 i = 1, nml 
IF(i.GE,j) GOTO 1 
T = Fa) 
m = F(i) 

F(i) = T 

1 k = nv2 

2 IF(k.GE.j) GOTO 3 
j = j-k 

k = k/2 
GOTO 2 

3 j = j + k 
DO 5 1 = 1, In 

le = 2**1 
lei = le/2 
U = (1.0,0.0) 

W = DCMPLX(DCOS(pi/lel),-DSIN(pi/lel)) 

DO 5 j = 1, lei 
DO 4 i = j,n,le 
ip = i + lei 
T = F(ip) * U 
F(ip) = F(i) - T 

4 F(i) = F(i) + T 

5 U = U • W 
RETURN 
END 
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SUBROUTINE FFT2D(DATA,TEMPDATA,F,ln,dir,n) 

c THIS S/R IS THE MAIN S/R AND PERFORMS A 2-D FFT ON AN 
c ARRAY OF SIDE LENGTH ’n’ WHERE Tn’ IS THE POWER OF 2. 
c ’dir’ IS EITHER +1 OR -1 WHETHER XFORMING OR INVERSE 

c XFORMING RESPECTIVELY. ’TEMPDATA’IS A WORKING ARRAY 
c FOR THE QUADRANT SWAPPING S/R. ’F IS THE 1-D ARRAY USED 
c BY THE 1-D FFT S/R. THIS S/R CALLS QUADSWAP, FFT, FOR BOTH 
c FORWARD AND INVERSE FFT AND CALLS NORMFFT AND 
c COMPLEXCONJ FOR INVERSE FFTs ONLY. THIS FFT NORMALIZES 
c BY DIVIDING BY n**2 WHEN THE INVERSE FFT IS PERFORMED. 

COMPLEX*16 DATA(n,n), F(n), TEMPDATA(n,n) 

INTEGER dir 

CALL Quadswap(DATA,TEMPDATA,n) 

IF (dir.EQ.-l) THEN 
CALL Complexconj(DATA,n) 

ENDEF 

DO 10 i = 1, n 
DO20j = l,n 
F(j) * DATA(ij) 

20 CONTINUE 
CALL FFT(F,ln,n) 

DO 30 j as 1, n 
DATA(ij) = FO) 

30 CONTINUE 
10 CONTINUE 
DO 40 j = 1, n 
DO 50 i = 1, n 
F(i) = DATA(ij) 

50 CONTINUE 
CALL FFT(F,ln,n) 
do 60 i = 1, n 
DATA(ij) = F(i) 

60 CONTINUE 
40 CONTINUE 

IF (dir.EQ.-l) THEN 
CALL Complexconj(DATA,n) 

CALL NormFFT(DATAn) 

ENDIF 

CALL Quadswap(DATA,TEMPDATA,n) 

RETURN 

END 
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SUBROUTINE Modulus(I,IS,IDATA,ISDATA,mod,icounter, 

+ nfiraines,nphoton,n) 

c THIS S/R DETERMINES THE MODULUS OF THE DEGRADED IMAGE 
c AND THE POINT SOURCE AND NORMALIZES THEM BY DIVIDING 
c BY THEIR RESPECTIVE DC VALUES. 


COMPLEX* 16 IDATA(-n/2:n/2-l,-n/2:n/2-l,0:4,-4:4), 
+ ISDATA(-n/2:n/2-l,-n/2:n/2-l), DC, DCS, 

+ I(-n/2:n/2-l,-ii/2:n/2-l), 

+ IS(-n/2:n/2-l,-n/2:n/2-l) 

REAL*8 mod(-ii/2:n/2-l,-n/2:n/2-l) 
k = 0 
ii2 = n/2 
n2ml = n2 -1 
DO 10 ii = -n2, n2ml 
DO 10 jj = -n2, n2ml 

IDATA(iijj,k,k) = IDATA(iijj,k,k) + ((I(iijj) * 

+ DCONJG(I(iiji))) - I(k,k)) 

ISDATA(ujj) = ISDATA(iiaj) + ((IS(iijj) * 

+ DCONJG(IS(iijj))) - IS(k,k)) 

10 CONTINUE 

IF (icounter.EQ.nframes) THEN 
DC = IDATA(k,k,k,k) 

DCS = ISDATA(k,k) 

DO 20 ii = -n2, n2ml 
DO 20 jj = -n2, n2ml 
IDATA(iijj,k,k) = IDATA(iijj,k,k)/DC 
ISDATA(iijj) = ISDATA(iiaj)/DCS 
ixiod(iijj) = nphoton * 

+ dsqrt(ABS(IDATA(ujj,k,k)/ISDATA(iioj))) 
IF (mod(ii jj).GT.dfloat(nphoton)) 

+ mod(iijj) = dfloat(nphoton) 

20 CONTINUE 
ENDIF 
RETURN 
END 
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SUBROUTINE NormFFT(DATA,n) 

c THIS S/R IS CALLED BY FTT2D FOR INVERSE FOURIER XFORMS 
c AND NORMALIZES THE XFORM BY DIVIDING BY n**2. 

COMPLEX* 16 DATA(n,n) 
nsqrd = n * n 
DO 10 i = 1, n 
DO 10 j = 1, n 
DATA(iJ) = DATA(ij)/nsqrd 
10 CONTINUE 
RETURN 
END 

SUBROUTINE Photons(DATA,nphoton,n) 

c THIS S/R NORMALIZES THE DATA ARRAY TO MAKE THE DC 
c VALUE EQUAL TO THE NUMBER OF PHOTONS IN THE OBJECT. 

COMPLEX* 16 DATA(n,n) 

REAL*8 photonum 
ic = n/2 + 1 

photonum = dfloat(nphoton)/DREAL(DATA(ic,ic)) 

DO 10 i = 1, n 
DO 10 j = 1, n 

DATA(ij) = photonum * DATA(ij) 

10 CONTINUE 
RETURN 
END 


SUBROUTINE Quadswap(DATA,TEMPDATA,n) 

c THIS S/R IS CALLED BY FFT2D AND SWAPS QUADRANTS OF DATA 

c ARRAY USING TEMPDATA AS A WORKING ARRAY. 

COMPLEX* 16 DATA(n,n), TEMPDATA(n,n) 
n2 = n/2 
DO 10 i = 1, n/2 
DO 10 j = 1, n/2 

TEMPDATA(ij) = DATA(i+n2j+n2) 

TEMPDATA(i+n2j) = DATA(i,j+n2) 

TEMPDATA(ij+n2) = DATA(i+n2j) 
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TEMPDATA(i+n2o+ii2) = DATA(ij) 
10 CONTINUE 
DO 20 i = 1, n 
DO20j = 1, n 

DATA(ij) = TEMPDATA(y) 

20 CONTINIJE 
RETURN 
END 


SUBROUTINE SNRcalc(KTsnr^nr,nyquist,n) 

c THIS S/R CALCULATES THE AVERAGE SNR AS A FUNCTION OF 
c RADIUS 

REAL*8 KTsiir(n,n), rsnr(n/2) 

INTEGER! 
n2 = ii/2 
n2pl = ii2 + 1 
DO 10 r = 1, nyquist 
nsnr = 0 
DO 20 i = 1, n 
DO20j - l,n 
X = float(j - (n2pl)) 
y = fioat(i - (n2pl)) 
radius = sqrt(x**2.0 + y**2.0) 

IF ((radius.GT.float(r-l))AND. 

+ (radius.LE.float(r))) THEN 
nsnr = nsnr + 1 
rsnr(r) = rsnr(r) + KTsnr(iJ) 

ENDIF 

20 CONTINUE 

rsnr(r) = rsnr(r)/nsnr 
10 CONTINUE 
RETURN 
END 
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SUBROUTINE Star(DATA,n) 

c THIS S/R CREATES A SIMULATED POINT SOURCE OR STAR 

COMPLEX*16 DATA(n,n) 
n2pl = i)/2 + 1 
DO 10 i = 1, n 
DO 10 j = l,n 
X = floatQ - n2pl) 
y = float(n2pl - i) 

IF ((x.EQ.0.0) AND.(y.EQ.0.0)) THEN 
DATA(id) = (1.0,0.0) 

ELSE 

DATA(ij) = (0.0,0.0) 

ENDIF 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE Writeffle(DATAl,data2,fflel,file2, 

+ nyquistyn) 

c THIS S/R WRITES THE DATA TO A FILE IN THE FORMAT REQUIRED 

c BY A PROGRAM WHICH PRODUCES A CONTOUR PLOT OF THE 
c IMAGE 

COMPLEX*16 DATAl(n,n) 

REAL*8 data2(n/2-l) 

REAL X, lambda, RO, pi 
INTEGER r 

CHARACIER*16 filel, ffle2 

COMMON A^ARS2/DIAM,OBSCUR,LAMBDA,RO,SECDIM 
PIXSCALE,TPFDIM,FILENAME 

OPEN (UNIT=30,FILE=filel,STATUS=’NEW’) 

OPEN (UNIT=40,FILE=file2,STATUS=’NEW’) 
pi = acos(-1.0E+00) 

DO 10 i = 1, n 
DO 10 j = l,n 
WRnE(30,*) DATAl(ij) 

10 CONTINUE 

DO 20 r * 1, nyquist 

X = r * (lamb^pixscale) * (180/pi) * 3600.0 
WRrrE(40,*) X, data2(r) 

20 CONTINUE 
RETURN 
END 
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APPENDK E. SPECIFIC KNOX-THOMPSON SUBROUTINES 


c THE FOLLOWING SUBROUTINES WERE GENERATED BY THE 
c THESIS AUTHOR AND ARE USED BY THE KNOX-THOMPSON 
c PROGRAM ONLY. 


SUBROUTINE UST 


SUBROUTINE Initialize(OBJDATA,GAUSSIAN,PUPIL,F, 

+ TEMPDATA,filel,file2,fwd4nv, 

+ lii,iiframes,nhoton,i^quist, 

+ ofEset,cent,n) 

c THIS S/R INITIALIZES SOME OF THE PROGRAM VARIABLES BY 
c QUEARYING THE USER FOR INPUT. IT ALSO CREATES THE 
c OBJECT DESIRED BY THE USER BY CALLING EITHER ASTEROID, 
c BISTAR, OR STAR S/R. 

INTEGER ofEset 
CHARACrER*16 filel, ffle2 
CHARACTER*! cent, object 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER INTEGER OFFSET (<= 4 AND <= #’ 
WRITE(*,*)TIXELS/RO):’ 

READ(*,*)of6set 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER INTEGER NUMBER OF SNAPSHOTS:’ 

READ(*,*)n£rames 

WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER INTEGER NUMBER OF OBJECT PHOTONS PER’ 
WRrTE(*,*)’SNAPSHOT:’ 

READ(*,*)nphoton 
WRTTEC*,*)’ ’ 

WRrTE(*,*)’ENTER INTEGER NYQUIST VALUE:’ 

READ(*,*)nyquist 
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10 WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER OBJECT TYPE: " A " FOR ASTEROID,’ 
Write(*,*)"' b " FOR BINARY STAR, OR " S " FOR STAR.’ 
READ(*,40) object 

IF ((object.EQ.’a’).OR.(object.EQ.’A’)) THEN 
CALL AST(OBJDATA,GAUSSIAN,PUPIL,F,TEMPDATA. 

+ fwd,iiiv4n,n) 

ELSEIF ((objcct.EQ.’b’).OR.(object.EQ.’B’)) THEN 
CALL Bistar(OBJDATA,n) 

ELSEIF ((object.EQ.’s’).OR.(object.EQ.’S’)) THEN 
CALL Star(OBJDATA,n) 

ELSE 

WRITE(*,*)’INCORRECT, REENTER’ 

GOTO 10 
ENDIF 

20 WRITE(*,*)’ ’ 

^'RITE(*,*)’WOULD YOU LIKE THE DATA CENTROIDED? (Y/N)’ 
READ(*,40) cent 

IF ((cent.NE.’y’)j\ND.(centNE.’Y’)AND. 

+ (cent.NE.’n’).AND.(cent.NE.’N’)) THEN 
WRITE(*,*)’INCORRECT, REENTER’ 

GOTO 20 
ENDIF 
WRITE(*,*)’ ’ 

WRITEC,*)’ENTER KT OUTPUT FILE NAME (16 CHAR’ 
WRITE(*,*)’MAX):’ 

READ(*,30) filel 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER KT SNR FILE NAME (16 CHAR’ 
WRriE(*,*)’MAX):’ 

READ(*,30) ffle2 
30 FORMAT(A16) 

40 FORMAT(Al) 

RETURN 

END 
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+ + 


SUBROUTINE KTrecon(KTPHASOR,I,IKT,KTsnr,offset, 
xvarTkt,xvaaikt,icounter, 
nframes 4 iyquist,n) 

c THIS S/R IS THE HEART OF THE KT RECONSTRUCTION PROGRAM, 

c IT DETERMINES THE VARIANCE-WEIGHTED, NOISE-BIAS- 
c CORRECHED (TIOSS-SPECTRUM OF THE DEGRADED IMAGE AND 
c RECONSTRUCTS THE PHASEORS FROM THE AVERAGE OF THIS 
c CROSS-SPECTRUM. THIS S/R IS CALLED ONCE FOR EACH 
c OBJECT/POINT SOURCE SNAP SHOT TO DETERMINE A RUNNING 
c AVERAGE CTIOSS-SPECTRUM AND WITH THE LAST SNAP SHOT, 
c THE PHASOR IS RECONSTRUCTED. 

COMPLEX*16 KTPHASOR(-n/2:ii/2-l,-n/2:n/2-l), CHEMP, 

+ IKT(-n/2:n/2-l,-n/2:n/2-l,-0:4,-4:4), 

+ I(-n/2:n/2-l,-n/2:n/2-l), PTEMPl, PTEMP2 
REAL*8 KTsnr(-n/2:ii/2-l,-n/2:n/2-l), sigmar, sigmai, 

+ xvarrlct(-n/2:ii/2-l,-ii/2:n/2-l,-0:4,-4:4), snr, 

+ xvarikt(-ii/2:ii/2-l,-n/2:n/2-l,-0:4,-4:4), sigma, 

INTEGER r, di, dj, offset 
k = 0 

KTPHASOR(k,k) = (1.0,0.0) 

DO 10 r = 1, nyquist 
DO 10 ii = 0, r 
IF (u.EQ.0) THEN 
lim = 0 
ELSE 
lim = -r 
ENDIF 

DO 10 jj = lim, r 

rad = sqrt(float(ii)**2.0 + float(jj)**2.0) 

IF ((rad.LE.float(r)).AND.(rad.GT.float(r-l))) THEN 
IF (icounter.EQ.n£rames) nsnr = 0 
DO 20 di = 0, offset 
DO 20 dj = -offset, offset 
drad = sqrt(float(di)**Z0 + float(dj)**2.0) 

IF ((drad.LE.float(r))AND. 

-I- (drad.LE.float(offset))) THEN 
idi = ii - di 
jdj = jj - dj 

radd = sqrt(float(idi)**2.0 + float(jdj)**2.0) 
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++ ++ ++ ++ ++ 


IF (raddXE.float(r-l)) THEN 
iHC = -ii 
jHC = -jj 

IF (I(idijdj).EQ.(0.0,0.0)) 

+ I(idijdj) = (1.0,0.0) 

IF (I(iijj).EQ.(0.0,0.0)) 

+ I(iijj) = (1.0,0.0) 

CTEMP = (I(idiodj) • DCONJG(I(iyi))) - 
+ DCX)NJG(I(di,dj)) 

IKT(idijdj,di,dj) = IKT(idydj,di,dj) + 

+ CTEMP 

xvarrkt(idijdj,di,dj) = 
xvarTkt(idijdj,di,dj) + 

(DREAL(CrEMP))**10 
xvarikt(idijdj,di,dj) = 
xvarikt(idydj,di,dj) + 

(DIMAG(CTEMP))**2.0 
IF (icounter.EQ.nframes) THEN 
nsnr = nsnr + 1 

sigmar = dabs((xvarrkt(idijdj,di,dj) - 
((DREAL(IKT(idijdj,di,dj)))**2.0)/ 
(^oat(nframes))/dfloat(nft^es • 1)) 
sigmai = dabs((xvarikt(idijdj,di,dj) • 
((DIMAG(IKT(idiJdj,di,dj)))**2.0)/ 
dfloat(iijErames))/dfloat(n&ames • 1)) 
sigma = dsqrt(sigmar + sigmai) 
snr = (((ABS(IKT(idiodj,di,dj)))/ 
dfloat(n£rames))/sigma) * 
dsqrt(dfloat(ii&^es)) 

KTsnr(iijj) = KTsnr(iijj) + snr 
PTEMPl = IKT(idijdj,di,dj)/ 

+ ABS(IKT(idijdj,di,dj)) 

PTEMP2 = KTPHASOR(idiJdj) 
KTPHASOR(ujj) = KTPHASOR(iijj) + 

+ ((snr**2.0) * (DCONJG(PTEMPl/PTEMP2))) 

ENDIF 
ENDIF 
ENDIF 

20 CONTINUE 
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IF (icounter.EQ.n£rames) THEN 
ICTsnr(iiji) = KTsnr(iijj)/nsnr 
KTsnr(iHCjHC) = KTsnr(iijj)/nsnr 
KTPHASOR(iiaj) = KTPHASOR(iiaj)/ 

+ ABS(KTPHASOR(iijj)) 

ICrPHASOR(iHCJHC) = DCONJG(KTPHASOR(iijj)) 
ENDIF 
ENDDF 

10 CONTINUE 
RETURN 
END 
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APPENDK F. SPECmC TRIPLECORREIATION SUBROUTINES 


c THE FOLLOWING SUBROUTINES WERE GENERATED BY THE 
c THESIS AUTHOR AND ARE USED BY THE TRIPLE-CORRELATION 
c PROGRAM ONLY. 


SUBROUTINE LIST 


SUBROUTINE BSrecon(BSPHASOR,I,IBS,IDBS,BSsnr,o£&et, 

-I- xvarrbs,xvaribs,icounter,nframes, 

+ i^quist,nphoton,n) 

c THIS S/R IS THE HEART OF THE BS RECONSTRUCTION PROGRAM, 

c rr DETERMINES THE VARIANCE-WEIGHTED, 
c NOISE-BIAS-CORRECTEDBISPECTRUM OFTHEDEGRADED IMAGE 

c AND RECONSTRUCTS THE PHASEORS FROM THE AVERAGE OF 
c THIS BISPECTRUM THIS S/R IS CALLED ONCE FOR EACH 
c OBJECT/POINT SOURCE SNAP SHOT TO DETERMINE A RUNNING 
c AVERAGE BISPECTRUM AND WITH THE LAST SNAP SHOT, THE 
c PHASOR IS RECONSTRUCTED. 

COMPLEX*16 BSPHASOR(-n/2:n/2-l,-ii/2;n/2-l), CTEMP, 

+ PTEMP2, IBS(-n/2:n/2-l,-n/2:n/2-l,0:4,-4:4), 

+ PTEMPl, I(-n/2:n/2-l,-n/2:n/2-l), ID^0:4,-4:4) 

REAL*8 BSsnr(-n/2:n/2-l,-n/2:n/2-l), sigmar, sigmai, 

+ snr, xvarrbs(-ii/2:n/2-l,-n/2:n/2-l,0:4,-4:4), 

+ sigma, xvaribs(-n/2:ii/2-l,-n/2:n/2-l,0:4,-4:4) 

INTEGER r, di, dj, offset 
k = 0 
loop = 1 

BSPHASOR(k,k) = (1.0,0.0) 

DO 10 r = 1, nyquist 
DO 10 ii = 0, r 
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+ + + + 


IF (ii.EQ.O) THEN 
lim s= 0 
ELSE 
lim = -r 
ENDIF 

DO 10 jj = lim, r 

rad = sqrt(float(ii)**2.0 + float(ij)**2.0) 

IF ((rad.LE.float(r))j\ND.(rad.GT.float(r-l))) THEN 
IF (icounter.EQ.iifiames) nsnr = 0 
DO 20 di = 0, of6set 
DO 20 dj = -offset, offset 
IF (loop.EQ.l) THEN 
IDBS(di,dj) = IDBS(di,dj) + I(di,dj) 

ENDIF 

drad = sqrt(float(di)**2.0 + float(dj)**2.0) 

IF ((drad.LE.float(r))AND. 

+ (drad.LE.float(offset))) THEN 
idi = ii - di 
jdj=jj-dj 

radd = sqrt(float(idi)**2.0 + float(jdj)**2.0) 

IF (radd.LE.float(r-l)) THEN 
iHC = -ii 
jHC = -jj 

IF (I(ididdj).EQ.(0.0,0.0)) 

+ I(idijdj) = (1.0,0.0) 

IF (I(iijj).EQ.(0.0,0.0)) 

+ I(iijj) = (1.0,0.0) 

CTEMP = I(idijdj) * 

+ DCONJG(I(iijj)) • 

+ I(di,dj) - (ABS(I(ididdj)))**2.0 - 

+ (ABS(I(iijj)))**2.0 - 

+ (ABS(I(di,dj)))**2.0 + (2.0 * nphoton) 

IBS(idiJdj,di,dj) = IBS(ididdj,di,dj) -f- 
+ CTEMP 

xvaiTbs(idijdj,di,dj) = 
xvarrbs(ididdj,di,dj) + 
(DREAL(CTEMP))**2.0 
xvaribs(idijdj,di,dj) = 
xvaribs(idijdj,di,dj) + 
(DIMAG(CrEMP))** 2.0 
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+ + + + + + 


IF (icounter.EQ.nframes) THEN 
nsnr = nsnr + 1 

sigmar = dabs((xvaiTbs(idiJdj,di,dj) - 
((DREA4IBS(idiodj,di,dj)))**2.0)/ 
dfloat(nframes))/dfloat(n£r^es - 1)) 
sigmai - dabs((xvaribs(idydj,di,dj) - 
((DIMAG(IBS(idijdj,di,dj)))**2.0)/ 
dfloat(nfirames))/dfloat(nfranies -1)) 
sigma - dsqrt(sigmar + sigmai) 
snr = (((ABS(IBS(idijdj,di,dj)))/ 
d£loat(nfimnes))/sigma) * 
dsqrt(dfloat(ii£rames)) 

BSsnr(iijj) = BSsnr(iyj) + snr 
PTEMPl = IBS(idiJdj,di,dj)/ 

+ ABS(IBS(idgdj,di,dj)) 

PTEMP2 = BSPHASOR(idi jdj) 

+ * IDBS(di,dj)/ABS(IDB^di,dj)) 

BSPHASOR(ujj) = BSPHASOR(iijj) + 

+ ((siir**2.0) * (DCONJG(PTEMPl/PTEMP2))) 

ENDIF 
ENDIF 
ENDIF 

20 CONTINUE 

IF (loop.EQ.l) loop as 2 
IF (icounter.EQ.nframes) THEN 
BSsnr(igj) = BSsnr(iijj)/nsnr 
BSsnr(iHCjHC) = BSsnr(iijj)/nsnr 
BSPHASOR(iijj) *= BSPHASOR(iyj)/ 

+ ABS(BSPHASOR(uJj)) 

BSPHASOR(iHCdHC) = DCONJG(BSPHASOR(iijj)) 
ENDIF 
ENDIF 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE Initialize(OBJDATA,GAUSSIAN,PUPIL,F, 

+ TEMPDATA,filel^e2,fwd,iiiv,lii, 

+ n&’ames,nphoton,nyqiiist,offset, 

+ cent,n) 

c THIS S/R INITIALIZES SOME OF THE PROGRAM VARIABLES BY 
c QUEARYING THE USER FOR INPUT. IT ALSO CREATES THE 
c OBJECT DESIRED BY THE USER BY CALLING EITHER ASTEROID, 
c BISTAR, OR STAR S/R. 

INTEGER offset 
CHARACTER*16 filel, file2 
CHARACTER*! cent, object 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER INTEGER OFFSET (<= 4):’ 

READ(*,*)offset 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER INTEGER NUMBER OF SNAPSHOTS:’ 

READ(*,*)n£rames 

WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER INTEGER NUMBER OF OBJECT PHOTONS PER’ 
WRITE(*,*)’SNAPSHOT:’ 

READ(*,*)nphoton 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER INTEGER NYQUIST VALUE:’ 

READ(*,*)nyquist 
10 WRITE(*,*)’ ’ 

WRnE(*,*)’ENTER OBJECT TYPE: " A " FOR ASTEROID,’ 
WRITEC*,*)"' B " FOR BINARY STAR, OR " S " FOR STAR.’ 

READ(*,40) object 

IF ((object.EQ.’a’).OR.(object.EQ.’A’)) THEN 

CALL ast(objdatagaussian,pupil,f,tempdata 

+ fwd,mv,ln,n) 

ELSEIF ((object.EQ.’b’).OR.(object.EQ.’B’)) THEN 
CALL Bistar(OBJDATAn) 

ELSEIF ((object.EQ.’s’).OR.(object.EQ.’S’)) THEN 
CALL Star(OBJDATA,n) 

ELSE 

WRrTE(*,*)’INCORRECT, REENTER’ 

GOTO 10 
ENDIF 
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20 WRITEC,*)’ ’ 

WRITE(*,*)’WOULD YOU LIKE THE DATA CENTROIDED? (Y/N)’ 
READ(*,40) cent 

IF ((ccnt.NE.y)AND.(ccnt.NE.*Y’)AND. 

+ (ccnt.NE.’n’)j\ND.(ccnt.NE.’N’)) THEN 
WRITE(*,*)’INCORRECT, REENTER’ 

GOTO 20 
ENDIF 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER TC OUTPUT FILE NAME (16’ 
WRITE(*,*)’(CHAR MAX):’ 

READ(*,30) filel 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER TC SNR FILE NAME (16 CHAR’ 
WRITE(*,*)’MAX):’ 

READ(*,30) file2 
30 FORMAT(A16) 

40 FORMAT(Al) 

RETURf4 

END 
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+ + 


APPENDIX G. ATMOSPHERIC DEGRADATION SUBROUTINES 


c THE FOLLOWING SUBROUTINES WERE PROVIDED BY CAPTAIN 
c CHUCK MATSON AND MS. IDA DRUNZER OF WL/ARCI AND 

c MODIFIED FOR USE WITH THE KNOX-THOMPSON AND 
c TRIPLE-CORRELATION IMAGE RECONSTRUCTION PROGRAMS, 
c THESE SUBROUTINES ARE REQUIRED FOR BOTH PROGRAMS. 


SUBROUTINE LIST 


SUBROUTINE PHSUB(SEIMG,F,TEMPDATA,OBJARR, 
icounter,nframes,nyquist, 
fwd,inv,b,n,mseed) 

c FUNCTION: Creates an image of a object located in space by filtering a 
c single snapshot of an object with a phase screen. 

c ORDER OF EVENTS TO OBTAIN IMAGE: 

c - Create an object 

c - Create Gaussian Array 

c - Expose array to what the atmosphere will do to the object 

c (Correlation Filter Function) 

c - Expose array to the stipulations of the telescope 

c - Autocorrelate the object (incoherent transfer function) 

c - Multiply the object by the phase screen to obtain the Image 

c DICTIONARY: 

c DIAM 
c ICOUNTER 
c 

c ISEED 
c 

c LAMBDA 
c MU 


- Telescope pupil diameter 

- Number of user chosen iterations of 

phase screens to process 

- Seed to be used for random number 

generator 

- Center wavelength 

- Mean 
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PDCSCALE - Corresponds to each pixel in (^cles/m 
RO - Variable characterizing the strength of atmospheric turbulence 

with respect to the optical system 
TPFDIM - Dimension for pupil of telescope 
PCOUNT - The actual number of photons per frame 
NFRAMES - The number of frames to process 
NPHOTON - The number of photons the user wants per frame 
NKL - The number of Karhunen-Loeve functions projected off and 

added back in the phase map. 

PARAMETER (IDIM = 64,LDIM = IDIM*IDIM) 

PARAMETER (IDIM2 = 64,LDIM2 = IDIM2*IDIM2) 

PARAMETER (IDIMl = .7071*IDIM2) 

PARAMETER (IWDIM = 5*IDIM/2,IWDIM2 = 5*IDIM2/2) 

PARAMETER (N1 = 11,N2 = (Nl/2) + 1, 

+ N3 = (N1 - 1)*(N1 + 2)/2) 

PARAMETER (N4 = (N1 - 1)*(N1 + 3)/4) 

PARAMETER (NKL = 20) 

CHARACTER*16 FILENAME 
CHARACTER*! TILT 

COMPLEX*16 PHAMAP(LDIM2),PS1(LDIM),IMAGE(LDIM), 

+ OBJARR(LDIM),PS2(LDIM2),SEIMG(LDIM) . 

REAL POIARRAY(LDIM),RVAR,MULTl,MULT2,LAMBDA, 

+ PIXSCALE,ZERO,SIZE,DIAM,TPFDIM,R0,INNERSCALE, 

+ CORR(LDIM2),SECDIM,OUTERSCALE,OBSCUR,RARRAY(LDIM) 

INTEGER PCOUNT 

DIMENSION ZKLMAP(IDIM2,IDIM2,NKL),LINPTR(N3),BVAL(N3) 

COMMON A^ARS/ZERO,RMAX,RMIN,CMINUSl,RCONSTl,FCONSTl,PI 
COMMON A^ARS2/DIAM,OBSCUR,LAMBDA,RO,SECDIM, 

+ PDCSCALE,TPFDIM,FILENAME 

c Initialize variables 

CMINUS1= (-l.E0,0.E0) 

ZERO = 0.E0 
MU = 0.E0 
SIGMA = 1.E0 


c 

c 

c 

c 

c 

c 

c 

c 
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PI = DACOS(-1.0D+00) 

JDIM = IDIM 
JDIM2 = IDIM2 
IC = IDIM/2+1 
INNERSCALE = .0001 
OUTERSCALE= 10000.E0 
MULTI = 5.92/INNERSCALE 
MULT2 = (2.0*PI)/OUTERSCALE 

c Only do certain subroutines the first time calling this subroutine (when wanting 
c several snapshots). 

11 CX)NTINUE 

IF((ICOUNTER.EQ.l)AND.(MSEED£Q.l)) THEN 
FILENAME = ’phvarld’ 

OPEN(10,FILE=FILENAME,STATUS=’OLD’,ERR=500) 

READ(10,*)DIAM,OBSCUR,LAMBDA,R0 

CLOSE(UNIT=10) 

c Call the subroutine which prompts the user for all names and variables that the 

c subroutine will need in order to run. 

CALL PARAM2(IDIM,NYQUIST,TILT,ISEED) 

c Generate a filter function to be used as a multiplier on a complex gaussian 

c array to be stored in array CORR. However, only call the filter function on 

c multiple runs if there are changing A^ue^ for the width of the object and RO. 

CALIFILT2(RARRAY,CORR,LDIM2,IDIM2,MULTl,MULT2,PIXSCALE, 
+ ZERO,RO) 

c Calculate the simulated size of the telescope (tpfdim). Compare it with the 
c dimension that is know (size). If size is less than tpfdim, there will be an error 

c due to the array being too small to store the data from the fourier transform, 

c and thus, data will be truncated. 

SIZE = (IDIM/2)-l 
IF(SIZE.LT.TPFDIM) THEN 
WRITE(*,*)’ ’ 

IF(SIZE.LE,(.5*TPFDIM)) THEN 
WRnE(*,*)’ERROR. IMAGE ARRAY SIZE MUST BE 
+ 512x512.’ 
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ELSE 

WRrrE(V)’ERROR. IMAGE ARRAY SIZE MUST BE 
+ 256x256.’ 


ENDEF 


WRITE(*,*)’ ’ 


GOTO 400 
ENDIF 


c Call BGSCREEN subroutine first time through phase screen. The BGSCREEN 
c subroutine manipulates the phase map to result in a phase map with proper 
c statistics. 


CALL BGSCR2(BVAL,ZKLMAP,LINPTR,IDIM2) 
c End the if statement if first time throu^ the subroutine, 
mseed = mseed + 1 


ENDIF 

c Get Gaussian numbers and then multiply these two numbers by the correlation 
c filter function to obtain the phase map (PHAMAP). 

DO 20 J = 1,LDIM2 

CALL GAUSS(MU,SIGMA,R1,R24>I,1SEED) 

PHAMAP(J) = DCMPLX(Rl*CORR(J),R2*CORR(J)) 

20 CONTINUE 

c Get phase screen which is in the frequency domain. Nrad perform the inverse 
c FFT on array. 

CALL FFT2D(PHAMAP,TEMPDATA,F,ln,inv,n) 

DO 25 J = 1, LDIM 

PHAMAP(J) = PHAMAP(J)*(DFLOAT(LDIM)) 

25 CONTINUE 

CALL PROJ(PHAMAP3VAL,ZKLMAP,LINPTR,iseed, 

+ NKL,R0,N1,N2,N3,N4,IDIM2^I,PDCSCALE,TILT) 
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c Calculate ae**(i*realization of phase screen). Take cosine and sine of eveiy 
c real element in the array and put into the phase screen array. 

DO 60 J = 1 JX)IM2 

PS2(J) = DCMPLX(DCOS(DREAL(PHAMAP(J))), 

+ DSIN(DREAL(PHAMAP(J)))) 

60 CONTINUE 

c Perform the Telescope Pupil Function (TPF) on the phase map to obtain the 
c coherent transfer function. First, call this routine to cut out chunk of phase 

c screen based on Telescope Pupil Function of D/LAMBDA to be stored in array 

c PSl. 

CALL TPF2A(PS14*S2,IDIM,IDIM2,TPFDIM,SECDIM) 

c Calculate the incoherent transfer function by multipfying phase screen by the 
c auto correlation filter function. Begin taking the Fourier Transform of the 
c phase screen. 

CALL FFT2D(PSl,TEMPDATA,F,ln,fwd,n) 

c Take the magnitude squared of the array and put the data into the real part of 
c the phase screen. 

DO 120 J = 1,LDIM 

PSl(J) = DCMPLX(DREAL(DCONJG(PSl(J))*PSl(J)),0.e0) 

120 CONTINUE 

CALL FFT2D(PSl,TEMPDATA,F,ln,inv,n) 

c Divide the phase screen by the value at dc in order to normalize the array by 
c that value (dcval). 

TEMP = DREAL(PS1(IC+(IC-1)*IDIM)) 

DO 125 J = 1,LDIM 
PS1(J) = PS1(J)/TEMP 
125 CONTINUE 
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c Incoherent transfer function is completed. Now get the image in the fourier 

c domain multiplying the object and the phase screen to get the image. 

DO 160 J = IJJDIM 
IMAGE(J) = (PSl(J)*OBJARR(J)) 

160 CONTINUE 

CALL FFr2D(IMAGE,TEMPDATA4?,ln4nv^) 

c Call the function, poisson, which will return as a floating point number in 
c integer value that is a random deviate drawn from a Poisson distribution of 
c mean equal to image value, using another real function, uniform, as a source • 

c of uniform random deviates. 

PCOUNT = ZERO 

DO 230 J = 1,LDIM 
RVAR = DREAL(IMAGE(J)) 

IF(RVAR.LE.0.0) THEN 
POIARRAY(J) = ZERO 

FT 

ITEMP = POISSON(RVAR,ISEED) 

POIARRAY(J) = DFLOAT(ITEMP) 

PCOUNT « PCOUNT + ITEMP 
ENDBF 

SEIMG(J) = DCMPLX(POIARRAY(J),ZERO) 

230 CONTINUE 

300 RETURN 
400 WRITE(V)” 

WRnE(*,*)’STOPPING PROGRAM DUE TO ERROR’ 

STOP 

c Error Statments 

500 CONTINUE 
WRITEC,*)’ ’ 

WRITE(V)’ERROR, THIS FILE DOES NOT EXIST. REENTER.’ 

WRITE(*,*)’ ’ 

GOTO 11 
800 CONTINUE 
END 
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c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


90 

95 


SUBROUTINE BGSCR2(bval,2klmap,linptr,nsid) 

This program does the things that result in a phase screen with proper statistics 


parameter(lu = 12,lu4 = 18) 
parameter(nr= 100^=20) 

PARAMETER (N1 = 11, N2 = (Nl/2) + 1, 
+ N3 = (N1 - 1)*(N1 + 1 ) 11 ) 

PARAMETER (N4 = (N1 - 1)*(N1 + 3)/4) 
real bval(n3),bvec(n2,n3),rmap(nr,n4),nsidl 
dimension zidmap(nsid,nsid,nbd) 
dimension ival(n3,3),ir(n3),linptr(n3) 
character* 16 foame,noname 
nsidl = .7071*(nsid -1) 
hiame = ’eigen.d’ 
noname = Idrad-d* 


Dictionary 


bval 

bvec 

ir 

ival 

linptr 

nl 

nkl 

nr 

nsid 

rmap 

zklmap 


contains eigenvalues of covariance matrix 

contains eigenvectors of covariance matrix 

pointers 

pointers 

pointers 

number of azimuthal orders to be used in the 2^mike functions 
number of Karhunen-l^oeve functions projected off and added 
back 

number of points in rmap (100 is more than enough) 
dimension of side of screen 

stored radial cut of Karhunen-Loeve functions, computed in 
michelin 

stored Karhunen-Loeve functions 


open(lu,file = fiiame,form=’formatted’,status=’old’) 
rewind(lu) 
read(lu,*) linptr 
rcad(lu,*) bval 
read(lu,*) ival 
do 95 i=l,n2 
do 90 j=l,n3 
read(lu,*) bvec(ij) 
continue 
continue 
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10 format(6E4.9) 

open(lu4,file»noname,form=’formatted*,status»’old’) 
rewind (lu4) 
do 80 i=l,nr 
do 85 j=l,n4 
rcad(lu4,*) rmap(ij) 

85 continue 
80 continue 
read(lu4,*) ir 

call kelp(zklmap,nnap,ir,ival,linptr^4i4, 

+ nsidynr,nkl) 

return 

end 


SUBROUTINE FADD(bmap,zklmap,iq,nsid,nkl,sum) 

c bmap is the phase screen 

complex* 16 bmap(nsid,nsid)^um 
dimension zklmap(nsid,nsid,n]d) 
data pi/3.14159265358979/ 

do 10 i = l,nsid 
do 20 j = l^id 

bmap(ij) « bmap(y) + sum*zklmap(y,iq) 

20 continue 
10 continue 

return 

end 


SUBROUTINE nLT2(RARRAY,C»RR,LDIM2,IDIM2,MULTl, 

+ MULT2,PIXSCALE®ROJRO) 

c Subroutine Function: This function generates an array representing the effects 
c of what the atmosphere will do to an object 

c when it passes through it. 

REAL RARRAY(LDIM2),CORR(LDIM2),MULTI,MULT2,RCONST, 

+ PIXSCALE,ZERO,RO,FCONST,R,ARRAYW 
INTEGER JOFFSET 
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c Calculate the radial average of an input image. 

c ARRAYS: CORR(LDIM2) - idim2 x idim2 array 
c RARRAY(IDIM2) > this is a small work array 

c that needs to be at least 

c IDIM2 in size 

c other variables: IDIM2 - x dimensionl 

c RCONST,RCONST2 - real constants 

c 

c load CORR with X (J) indicies squared 

FCONST = .1517 
DO 10 J=1,IDIM2 
JOFFSET=(J-l)*IDIM2 
R=FLOAT( J-(IDIM2/2+l) ) 

RCONST=R*R 
DO 15 I=1,IDIM2 
CORR(I+JOFFSET) = RCONST 
15 CONTINIJE 
10 CONTINUE 

c Add to LI, Y (I) indicies squared (one column of Y indicies stored in L2) 

DO 20 I=1,IDIM2 
R=FLOAT( I.(IDIM2/2+l)) 

RARRAY(I) = R*R 
20 CONTINUE 
DO 30 J=1,IDIM2 
JOFFSET=(J-l)*IDIM2 
DO 40 I = 1,IDIM2 

CORR(JOFFSET + I) = RARRAY(I) + CORR(JOFFSET + I) 

40 CONTINUE 
30 CONTINUE 

c Move the scaling array to another array in order to do calculations for the filter 

c function 

DO 50 I = 1,LDIM2 
RARRAY(I) = CORR(I) 

50 CONTINUE 
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c Multiply the scaling array by the correlation filter function 


RCX)NST = (-11^/12.E0) 

DO 601 = 

RARRAY(I) = (RARRAY(I)+(MULT2*MlJLT2))**RCONST 
60 CONTINUE 

ARRAYW = PIXSCALE*IDIM2 
RCONST = (R0/ARRAYW)**(-5.E0/6.E0) 

DO 701 = 1,LDIM2 

RARRAY(I) = (RARRAY(I)*RCONST)*FCONST 
70 CONTINUE 

c Set the middle point in array to zero to normalize the array by that value 

RARRAY(IDIM2*(IDIM2/2)+(IDIM2/2+l)) = ZERO 

c Calculate second part of correlation filter function and multiply it the first 
c part to complete the filter array. 

RCONST = -2.0*(MULT1*MULT1) 

DO 90 1= 1,LDIM2 

CORR(I) = (EXP(CORR(I)/RCONST))*RARRAY(I) 

90 CONTINUE 

900 CONTINUE 
RETURN 
END 


SUBROUTINE GAUSS(MU,SIGMA,RNUMl,RI»UM2,PI,iseed) 
c PURPOSE: 

c GET GAUSSIAN DISTRIBUTED RANDOM NUMBER WITH MEAN MU 

c AND STANDARD DEVIATION SIGMA 

REAL MU,SIGMA,Yl,Y2,RNUMl,RNUM2,PI,TWOPI 

TWOPI = 2.eO*PI 

Yl= uniform(iseed) 

Y2= uniform(iseed) 
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c Calculate first random number 

RNUMl = MU + SIGMA*SQRT(-2.0*ALOG(Y1))*COS(TWOPI*Y2) 

c Calcidate the second random number 

RNUM2 = MU + SIGMA*SQRT(-ZO*AIjOG(Yl))*SIN(TWOPrY2) 

RETURN 

END 


SUBROUTINE INPROD(bmap;eklmap,iq,nsid,nkl^p) 

complex* 16 bmap(nsid,nsid) 

complex* 16 zap 

dimension zklmap(nsid,nsid,nkl) 

zap = (0.0,0.0) 
area = 0.0 
do 10 i = l,nsid 
do 20 j = l,nsid 

zap = zap + bmap(ij)*zklmap(ij,iq) 
area *= area + ( zklmap(i j,iq) )**2 
20 continue 
10 continue 

zap = zap/area 

return 

end 
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SUBROUTINE KELP(zklmap^ap,ir,ival,lmptr, 

c zklmap is the map of the k-1 funcdons, sigh 

dimension rmap(nr,n4),ir(n3) 
dimension zkhnap(nsi^nsid,nkl) 
dimension ival(n3,3)4mptr(n3) 
data pi/3.14159265358979/ 
icent = nsid/2 + 1 
dx = 1.0/(nsid/2 -1) 
dr = 1.0/(nr -1) 
do 15 iq = 1^ 
m * ival( linptr(iq),l ) 
iod = ival( linptr(iq),2 ) 
ipq = ir(iq) 
do 10 i = l,nsid 
X = (i - icent)*dx 

do 20 j s l^isid 

zkhnap(ij,iq) = 0.0 

y = (j - icent)*dx 
r = sqrt(x**2 + y**2) 
th = 0.0 
if( r.le.1.0) then 
if( r.gt.0.0) then 
th = atan2(y^) 
if(th.le.0.0) th = 2.0*pi + th 
else 

th = 0.0 
endif 

if(r.lt.0.99999) then 
kl = int(r/dr) + 1 
k2 = kl + 1 
rl = (kl - l)*dr 
r2 = (k2 - l)*dr 
coef 1 = (r2 - r)/dr 
coe£2 = (r - riydr 

top = coefl*rmap(kl,ipq) + coef2*rmap(k2,ipq) 
else 

top = rmap(nr,ipq) 
endif 

if(iod.eq.l) zz = cos(m*th) 
if(iod.eq.2) zz = sin(m*th) 
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zklmap(y,iq) = 22 *top 
endif 
continue 
continue 
continue 
return 
end 


SUBROUTINE PARAM2(IDIM^QUIST,TILT,ISEED) 

The function of this subroutine is to ask the user various questions about the 
subroutines he or she wishes to use and what values he or she wants to assign 
to the variables in the program. 

REAL DIAM,LAMBDA,R0,PIXELNUM,OBSCUR4*IXSCALE, 

+ secdihtpfdim 

INTEGER IDIM 
CHARACTER* 16 FILENAME 
CHARACTER*! ANSWER,FLAG,TILT 
COMMONA^ARS2/DIAM,OBSCUR,LAMBDA,RO,SECDIM, 

+ PIXSCALE,TPFDIMJTLENAME 

WR1TE(*,*)’ ’ 

WRITE(*,*)’THE DEFAULT VALUES FOR THIS PROGRAM ARE AS 
+ FOLLOWS: ’ 

WRITE(*,*)’ ’ 

WRrTE(*,*)’ TELESCOPE DIAMETER: ’,DIAM 

WRITE(*,*)’ CENTER WAVELENGTH: LAMBDA 

WRrTE(*,*)’ RO: ’,R0 

WRrTE(*,*)’ TELESCOPE SECONDARY DIAMETER:’,OBSCUR 

CONTINUE 
WRITEC*,*)’ ’ 

WRrTE(*,*)’USE THE DEFAULT VALUES? (Y/N)’ 

READ(*,14) FLAG 

IF(((FLAG.NE.’Y’)J^,(FLAG.NE.’y’))j\ND. 

+ ((FLAG.NE.’N’).AND.(FLAG.NE.’n’))) THEN 
WRrTE(*,*)’ ’ 

WRITE(*,*)TOUR ANSWER HAS TO BE EITHER Y OR N. 

+ REENTER’ 

FLAG = ” 

GOTO 20 
ENDIF 
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30 CONTINUE 

IF((FLAG.EQ.’N’).OR.(FLAG.EQ.’n’))THEN 
WRITEC*,*)’ ’ 

WRITE(V)’ A: ALL VARIABLES’ 

WRITEC*,*)’ D: TELESCOPE DIAMETER’ 

WRITE(*,*)’ W: CENTER WAVELENGTH (LAMBDA)’ 
WRnE(*,*)’ R: RO’ 

WRITE(*,*)’ S: SECONDARY DIAMETER’ 

WRITE(*,*)’ ’ 

WRITE(*,*)’WHICH VALUE WOULD YOU LIKE TO CHANGE?’ 
READ(*,14) ANSWER 

IF(((ANSWER.NE.’A’)AND.(ANSWER.NR’a’))j\ND. 

+ ((ANSWER.NE.’D’)J^ND.(ANSWER.NE.’d’))j\ND. 

+ ((ANSWERJ^.’L’)AND.(ANSWER.NE.’1’))AND. 

+ ((ANSWER.NE.’W’).AND.(ANSWER.NE.’w’))AND. 

+ ((ANSWEILNE.’R’)AND.(ANSWER.NE.’r’))AND. 

+ ((ANSWER.NE.*S’)AND.(ANSWER.NE.’s’)))THEN 
WRITE(*,*)’ ’ 

WRITER,*)’ANSWER IS NOT A CHOICE. REENTER’ 

GOTO 30 
ENDIF 
WRnE(*,*)’ ’ 

WRITEC*,*)*!™ CURRENT VALUES ARE: ’ 

WRITE(* •)’ ’ 

WRITE(*,*)’ TELESCOPE DIAMETER: ’,DIAM 

WRnE(*,*)’ CENTER WAVELENGTH: ’,LAMBDA 

WRnE(*,*)’ RO: ’,R0 

WRITE(*,*)’ SECONDARY DIAMETER ’,OBSCUR 

WRnE(*,*)’ ’ 

WRITE(*,*)’ ’ 

IF((ANSWER.EQ.’A’).OR.(ANSWER.EQ.’a’))THEN 
WRnE(*,*)’INPUT THE FOLLOWING’ 

WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER NEW TELESCOPE DIAMETER:’ 
READ(*,100) DIAM 
WRrrE(*,*)’ ’ 

WRITE(*,*)’ENTER NEW CENTER WAVELENGTH:’ 
READ(*,100) LAMBDA 
WRITEC,*)’ ’ 

WRITE(*,*)’ENTER NEW RO:’ 

READ(*,100) RO 
WRITEC*,*)’ ’ 

WRITER,•)’ENTER NEW TELESCOPE SECONDARY DIAM:’ 
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50 


READ(*,100) OBSCUR 

ELSEIF((ANSWER.EQ.’D’).OR.(ANSWER.EQ.’d’))THEN 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER NEW TELESCOPE DIAMETER:’ 
READ(*,100)DIAM 

ELSEIF((ANSWER.EQ.’W’).OR.(ANSWER.EQ.’w’))THEN 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER NEW CENTER WAVELENGTH:’ 
READ(*,100) LAMBDA 

ELSEIF((ANSWER.EQ.’R’).OR.(ANSWER.EQ.’r’))THEN 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER THE NEW RO:’ 

READ(*,100)R0 

ELSEIF((ANSWER.EQ.’S’).OR.(ANSWER.EQ.’s’))THEN 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER THE NEW SECONDARY DIAMETER:’ 
READ(*,100) OBSCUR 
ENDEF 
CONTINUE 


WRrrE(*,*)’ ’ 

WRITER,•)’WOULD YOU LIKE TO CHANGE ANY OTHER VALUE? 
+ (Y/N)’ 

READ(*,14)ANSWER 

IF(((ANSWER.NE.’Y’).AND.(ANSWER.NE.’y’)).AND. 

+ ((ANSWER.NE.’N’).AND.(ANSWER.NE.’n’)))THEN 


WRITE(*,*)’ ’ 

WRnE(*,*)’YOUR ANSWER HAS TO BE EITHER Y OR N. 
+ REENTER’ 

GOTO 50 
ENDIF 

IF((ANSWER.EQ.’Y’).OR.(ANSWER.EQ.’y’))THEN 
GOTO 30 
ENDIF 

OPEN(UNIT= 10,FILE=FILENAME,FORM=’FORMATTED’, 

+ STATUS=’UNKNOWN’) 

REWIND(IO) 

WRITE(10,*)DIAM,OBSCUR,LAMBDA,R0 
CLOSE(UNIT=10) 

ENDIF 


WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER THE NUMBER OF PIXELS PER RO:’ 
READ(*,*)PDCELNUM 
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WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER SEED VALUE:* 

READ(*,*)ISEED 
PDCSCALE = RO/PDCELNUM 
TPFDIM » nmt(DIAM*PDCELNUM/RO) 

SECDIM = iimt(OBSCUR*PDCELNUM/RO) 

60 CONTINUE 
65 CONTINUE 
WRITE(*,*)’ ’ 

WRITE(*,*)’IX) YOU WISH TO INCLUDE X/Y TILT? (Y/N)’ 
READ(*,14)TILT 

IF(((TiTJ^.’Y’)j\ND.(TILT.NE.y))AND. 

+ ((TILT.NE.’N’)j\ND.(TILT.NE.’n’)))THEN 
WRITE(*,*)’ ’ 

WRITE(*,*)*YOUR ANSWER MUST BE EITHER Y OR N. 
+ REENTER’ 

GOTO 65 
ENDIF 

c Check the value for nyquist to assure it is equal to tpfdim 
IF(NYQUIST.NE.TPFDIM) GOTO 999 
RETURN 

c Format Statements 

14 FORMAT(Al) 

100 FORMAT(E11.4) 
no FORMAT(I6) 

c Error statements 

999 WRITE(V)” 

WRITE(*,*)’ERROR. VALUE FOR NYQUIST MUST EQUAL 
+ TO:’,TPFDIM 

+ ,’ NOT ’,nyquist, ’.’ 

WRITE(*,*)’ ’ 

WRirE(*,*)TROGRAM STOPPING. CHANGE NYQUIST.’ 

STOP 

END 
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SUBROUTINE PROJ(bmap,bv al,zklnia p,linptr,iseed, 

+ nkl,r0,nl,n2,n3,n4,idiin2,pi,pixs(^e,tilt) 

c The function of this subroutine is to project off the low-order Karhunen-Loeve 
c functions and then adds them back with the proper strength. 

c nl is the number of azimuthal orders, bmap is the phase screen 

complex* 16 bmap(idim2,idim2),zap,sum 
dimension bval(n3),linptr(n3) 
dimension zklmap(idim2,idii]^nkl) 
real x,y,yi,pi,twopi,pixscale,drO,rO 
integer izem,idi^ 
character tilt 


izem = 5 

twopi = 2.e0*pi 

drO = (pixscale*idim2)/r0 

DO 10 i = l,izem 

c A random size is generated for each karhunen function 

X = sqrt(bval(linptr(i) ) ) 
y = O.eO 
yi = 0.e0 
mu = 0.e0 
sigma = X 

call gauss(mu,sigma,y,yi,pi,iseed) 

sum = cmplx(Y,YI) 

sum = sum*(dr0* *.83333333) 

call inprod(bmap,zklmap,i,idim2,nkl,zap) 

c The function is projected out 

if(tilt.eq.’N’.or.tilt.eq.’n’) then 
if(i.le.2)sum = 0.0 
endif 

sum = sum - zap 

call fadd(bmap,zklmap,i,idim2,nkl,sum) 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE TPF2A(PS1,PS2,IDIM,IDIM2,TPFDIM,SECDIM) 

c Function: The function of this subroutine is to cut out a portion of array with 
c the dimension of TPFDIM and put the data into an array of zeros. 

COMPLEX*16 PS1(IDIM,IDIM),PS2(IDIM2,IDIM2) 

REAL TPFDIM,SECDIM 
INTEGER IC,IC2 

IC = IDIM/2 + 1 
IC2 = IDIM2/2 + 1 

c Zero out psl to ensure that no reminents of previous phase screen is left over. 

DO 10 J = 1,IDIM 
DO 10 K = 1,IDIM 
PS1(J,K) = (0.0,0.0) 

10 CONTINUE 

c Divide the telescope pupil dimension 2 in order to obtain the radius instead 
c of the diameter of the pupil. 

TPFDIM2 * TPFDIM/2 
SECDIM2 * SECDIM/2 

IF(SECDIM.EQ.0.0)PS1(IC,IC) = PS2(IC2,IC2) 

c Now insert the phase screen into the larger array. The smaller array is the data 

c that this telescope is able to see due to the inhibition of its size. 

DO 20 I = 0, TPFDIM2 
DO 20 J * 0, TPFDIM2 
IF((r*2 + J**2).LE.TPFDIM2**2) THEN 
IF((I**2 + J**2).GT.SECDIM2**2) THEN 
PS1(IC+I,IC+J) = PS2(IC2+I,IC2+J) 

PS1(IC+I,IC-J) = PS2(IC2+I,IC2-J) 

PS1(1C-I,IC+J) = PS2(1C2-I,IC2+J) 

PS1(IC-I,IC-J) * PS2(IC2-I,IC2-J) 

ENDIF 

ENDIF 

20 CONTINUE 
RETURN 
END 


158 




FUNCTION GAMMLN(xx) 

real cof(6),stp,half^one4pf,x,tmp^r 

data coUtp^6.18009173d0,-86.50532033d0, 

+ 24.01409822dO,-1.231739516dO,.120858003d-2, 
+ -J36382d-5,2J0662827465d0/ 
data half,one,^^.Sd0,l>0d0,5.5d0/ 


X = XX - one 

tmp = X + 

tmp = (x + half)*log(tmp) - tmp 
ser = one 

do lOj = 1,6 
X = X + one 
ser = ser + cof(j)/x 
10 continue 

gammln = tmp + log(stp*ser) 

return 

end 


REAL FUNCTION POISSON(pmean, iseed) 
c mean of poisson distribution 
real pmean 

c Seed for random number generator 
integer iseed 

c Summary of purpose 

c Generate a random number with a poisson distributionof mean pmean (poisson 

c deviate) 

c Author 

c numerical recipes, p207, routine poidev.for 
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c Modifications 
c 13 nov, 1987: pat fitch 

c 1988: erik m Johansson - fixed integer conversion problem with maxint 

c Routines called 

c uniform 
c gammln 

c (c) Copyright 1987 the Regents of the University of California. All rights 
c reserv^. 

c This software is a result of work performed at Lawrence Livermore National 
c Laboratory. The United States Government retains certain rights therein. 

real pi 

parameter (pi=3.141592654) 

c maxint is the largest real number which can be converted to an integer without 
c resulting in an arithmetic error (32 bits) 

real maxint 

parameter (maxint = .214748352el0) 

real pexpmean, oldmean, t, em, sq, alxm^ y, gammln, 

+ uiiiform 
integer times 
data oldmean Al./ 
if(pmean.lt. 12.0)then 

c Use direct computation method 

if(pmean.ne.oldmean) then 

c New mean, calculate the exponential 

oldmean = pmean 
pexpmean= exp(-pmean) 
endif 
em = -1.0 
t = 1.0 

2 em = em + 1.0 
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t = t * iiniform(iseed) 
if(tg:t.pexpinean) go to 2 
else 

if(pmeanjie.oldmean) then 
oldmean = pmean 
sq = sqrt(2.0*pmean) 
ato= alog(pmean) 

c Natural log of gamma function = gammln 

pexpmean= pmean*alxm - gammln(pmean+l.) 
endif 


times = 0 

1 y = tan( pi* uniform(iseed)) 
times = times + 1 
if (times .ge. 1000) then 

write(*,*)’ERROR: STUCK IN lX)OP IN POISSON’ 
stop 
endif 

em = sq*y+pmean 

if((em .It. 0.) .or. (em .gt. maxint)) go to 1 
em = int(em) 

t « 0.9*(l.+y**2)*exp(em*alxm-gammln(em+l.)- 
+ pexpmean) 
if(uniform(iseed).gtt) go to 1 
endif 

poisson = em 

return 

end 
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REAL FUNCTION UNIFORM(seed) 
integer*4 seed 
c Summary of purpose 

c implements a multiplicative linear congruential generator generates random 
c numbers uniformly distributed on the open interval 0.0 to 1.0 (0 and 1 are NOT 
c included) 

c Author 

c from "Random Number Generators: Good Ones Are Hard to Find," Stephen 
c K. Park and Keith W. Miller, Co mmuni cations of the ACM, October 1988, Vol 
c 31, No 10, pp 1192 - 1201 routine is listed on p 1195 

c modifications 

c 4/11/89 erik m Johansson - modified code to use less computational steps, 
c Equations used are from the article "Efficient and Portable Combined Random 

c Number Generators," Pierre L’Ecuyer, Communications of the ACM, June 
c 1988, Vol 31, No 6, at the top of p745. 

c (c) Copyright 1989 the Regents of the University of California. All rights 
c reserved. 

c This software is a resxilt of work performed at Lawrence Livermore National 
c Laboratory. The United States Government retains certain rights therein. 

integer*4 a, m, q, r 

parameter (a = 16807) 
parameter (m = 2147483647) 
parameter (q = 127773) 
parameter (r = 2836) 

real h 

parameter (h = 1. / 2147483647.) 
integer*4 k 
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c The function 


k = seed / q 

seed = a * (seed - k*q)-k*r 

if (seed .It. 0) seed = seed + m 

uniform = seed * h 

return 

end 
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APPENDK H. PHASE ERROR AND LOW PASS FILTER PROGRAM 


c THIS PROGRAM WAS DEVELOPED BY THE THESIS AUTHOR AND 
c DEVELOPS A LOW PASS FILTER IN THE FREQUENCY 
c DOMAIN TO FILTER THE RECONSTRUCTED OBJECT DATA FILE 
c FOR BOTH THE KNOX-THOMPSON AND TRIPLE-CORRELATION 
c METHODS. ADDITIONALLY, THIS PROGRAM DETERMINES THE 
c AZIMUTHAL RMS PHASE ERROR FOR THE INPUT FOURIER 
c SPECTRA. 

c THE FOLLOWING SUBROUTINES ARE REQUIRED FROM 
c UNIVERSAL SUBROUTINES IN APPENDDC D: 

c Complexconj 

c FFT 

c FFT2D 

c NormFFT 

c Quadswap 


c AUTHOR: 
c COMPL. DATE: 
c REASON: 
c 

c GOAL: 
c 
c 
c 


LT JAMES M. LACKEMACHER 
26 OCTOBER 1990 

COMPLETE REQUIREMENTS FOR A MASTERS 
DEGREE IN PHYSICS 
SIMULATE OBJECT, DEGRADE OBJECT, 
RECONSTRUCT OBJECT USING KNOX-THOMPSON 
AND TRIPLE-CORRELATION METHODS FILTER 
THE DATA THE COMPARE THE TWO METHODS. 


PROGRAM FREQFILTER 

MAIN PROGRAM COMPLEX VARIABLE LIST 

c BSDATA n X n DIM ARRAY REPRESENTING THE INPUT 

c TRIPLE-CORRELATION RECONSTRUCTED OBJECT IN 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BSLP 

F 

KTDATA 

KTLP 

TEMPDATA 

TRUDATA 


FREQUENCY SPACE 

n X n DIM ARRAY REPRESENTING THE LOW PASS 
FILTERED TRIPLE-CORRELATION RECONSTRUCTED 
OBJECT IN FREQUENCY SPACE 
n DIM ARRAY USED IN THE FOURIER TRANSFORM 
n X n DIM ARRAY REPRESENTING THE INPUT 
KNOX-THOMPSON RECONSTRUCTED OBJECT IN 
FREQUENCY SPACE 

n X n DIM ARRAY REPRESENTING THE LOW PASS 
FILTERED KNOX-THOMPSON RECONSTRUCTED 
OBJECT IN FREQUENCY SPACE 
n X n DIM ARRAY THAT IS USED AS A TEMPORARY 
ARRAY IN THE FOURIER TRANSFORM 
n X n DIM ARRAY REPRESENTING THE TRUTH DATA 


c 


MAIN PROGRAM REAL VARIABLE LIST 


r 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


bserror 

bslpmod 

bsr 

bssnr 

bstninc 

kterror 

ktlpmod 

ktr 

ktsnr 


n X n DIM ARRAY THAT REPRESENTS THE SQUARE 
PHASE ERROR FOR THE TRIPLE-CORRELATION 
RECONSTRUCTED IMAGE 

n X n DIM ARRAY THAT REPRESENTS THE MODULUS 
OF THE LOW PASS FILTERED TRIPLE-CORRELATION 
RECONSTRUCTED IMAGE 

n/2 DIM ARRAY THAT REPRESENTS THE RADIAL 
FREQ VALUE IN INVERSE ARCSECONDS FOR THE 
TRIPLE-CORRELATION IMAGE 
n/2 DIM ARRAY THAT REPRESENTS THE INPUT SNR 
VALUES OF TRIPLE-CORRELATION IMAGE 
TRIPLE-CORRELATION RADIALTRUNCATION VALUE 
n X n DIM ARRAY THAT REPRESENTS THE SQUARE 
PHASE ERROR FOR THE KNOX-THOMPSON 
RECONSTRUCTED IMAGE 

n X n DIM ARRAY THAT REPRESENTS THE MODULUS 
OF THE LOW PASS FILTERED KNOX-THOMPSON 
RECONSTRUCTED IMAGE 

n/2 DIM ARRAY THAT REPRESENTS THE RADIAL 
FREQ VALUE IN INVERSE ARCSECONDS FOR THE 
KNOX-THOMPSON IMAGE 

n/2 DIM ARRAY THAT REPRESENTS THE INPUT SNR 
VALUES OF KNOX-THOMPSON IMAGE 
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c kttnmc 
c lambda 
c rBSerror 
c 
c 

c rKTeixor 

c 

c 

c RO 


KNOX-THOMPSON RADIAL TRUNCATION VALUE 
WAVELENGTH OF QUASI-MONOCHROMATIC UGHT 
n/2 DIM ARRAY THAT REPRESENTS THE AZIMUTHAL 
RMS PHASE ERROR OF THE TRIPLE-CORRELATION 
RECONSTRUCTED IMAGE 

n/2 DIM ARRAY THAT REPRESENTS THE AZIMUTHAL 
RMS PHASE ERROR OF THE KNOX-THOMPSON 
RECONSTRUCTED IMAGE 
COHERENCE LENGTH 


c MAIN PROGRAM INTEGER VARIABLE UST 


c fwd 
c iktcount 
c inv 
c itccount 
c In 
c n 
c nyquist 
c 
c 

c numpix 


VALUE OF 1 FOR FORWARD FFT 

NUMBER OF KT PIXELS WITH SNR GREATER THAN 1.0 

VALUE OF -1 FOR INVERSE FFT 

NUMBER OF TC PIXELS WITH SNR GREATER THAN 1.0 

2'"In FOR USE WITH FFT SUBROUTINE 

DIMENSION OF ONE SIDE OF 2-DIM ARRAY 

EQUAL TO THE TELESCOPE PUPIL FUNCTION 

DERIVED FROM THE FOLLOWING FORMULA: 

nyquist = (telescope diameter x number of pixels/r0)/r0 

THE NUMBER OF PIXELS PER RO 


c MAIN PROGRAM INTEGER VARIABLE LIST 


c 

end 

CHAR INPUT TO DETERMINE WHETHER TO STOP THE 

c 


PROGRAM 

c 

mean 

CHAR INPUT TO DETERMINE WHETHER TO FIND THE 

c 


AMSPE 

c 

rite 

CHAR INPUT TO DETERMINE WHETHER TO WRITE 

c 


DATA TO FILE 


f 


166 





MAIN PROGRAM 


PARAMETER(n=644n=6,fwd= l,inv=-l) 


COMPLEX* 16 KTDATA(n,n), BSDATA(n,n), TEMPDATA(n,n), 
+ F(n), KTLP(n,n), BSLP(n,n), TRUDATA(n,n) 

REAL*8 ktlpmod(n,n), bslpmcxl(n,n), kterror(n,n), 

+ bseiTor(n,n), rCTerror(n/2), rBSeiTor(ii/2), 

+ ktsnr(ii/2), bssnr(n/2) 

REAL ktr(n/2), bsr(i)/2), kttrunc, bstrunc, lambda 
CHARACTER*! mean, end, rite 


c INPUT THE DATA 


CALL Readffle(KTDATA,BSDATA,n) 

WRITE(*,*)’ENTER NYQUIST VALUE:’ 

READ(*,*)nyquist 
WRin^*,*)’ ’ 

WRITE(*,*)’ENTER THE NUMBER OF PDCELS PER RO:’ 
READ(*,*) numpix 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER RO VALUE IN METERS:’ 

READ(*,*) RO 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER WAVELENGTH VALUE IN METERS:’ 
READ(*,*) lambda 


WRITE(*,*)’ ’ 

CALL Readsnr(ktsnr,bssnr,ktr,bsr,nyquist,n) 


c LOW PASS FILTER BASED ON RADIAL TRUNCATION VALUE 
c DETERMINED FROM SNR DATA 


CALL Tnincval(ktsnr,ktr,kttrunc,nyquist,iktcount,n) 
WRrrE(*,*)’KT RADIAL TRUNCATION VALUE IS:’ 
WRITE(*,*) kttnmc,’l/ARCSEC.’ 

WRITE(*,*)’VALUES GREATER THAN THIS RADIUS’ 
WRITE^,*)’WILL BE TRUNCATED.’ 

WRITEC*,*)’ ’ 
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CALL Truncval(bssnr,bsr,bstrunc,nyquist,itccount,n) 

WRrrE(*,*)TC RADIAL TRUNCATION VALUE IS:’ 

WRITE(*,*) bstrunc,’l/ARCSEC’ 

WRITE(*,*)’VALUES GREATER THAN THIS RADIUS’ 
WRITE(*,*)’WILL BE TRUNCATED.’ 

WRITE(*,*)’ ’ 

c TRUNCATE THE DATA 

CALL Tnincate(KTDATA,KTLP,kttruiic,lainbda,RO,numpix,n) 
CALLTruncate(BSDATA,BSLP,bstrunc,lambda,RO,numpix,n) 

c INVERSE FFT THE DATA AFTER FILTERING 

CALL FFT2D(KTLP,TEMPDATA,F,ln,mv,n) 

CALL FFT2D(BSLP,TEMPDATA,F,ln,mv,n) 

c DETERMINE THE MODULUS OF THE DATA IN IMAGE SPACE 

CALL Modulus(ktlpmod,KTLP,n) 

CALL Modulus(bslpmo(l,BSLP,n) 

c NORMALIZE THE MODULI FOR PLOTTING 

CALL Normalize(ktlpmod,n) 

CALL Nomialize(bslpmod,n) 

c WRITE THE MODULUS TO A FILE FOR PLOTTING IF DESIRED 

20 WRITER,•)’WRnE MODULUS TO A FILE? (Y/N)’ 

READ(*,*)rite 
WRITE(*,*)’ ’ 

IF ((rite.NE.’Y’)j\ND.(rite.NE.’y’).AND. 

+ (rite.NE.’N’).AND,(rite.NE.’n’)) THEN 
WRrTE(*,*)’ERROR, REENTER.’ 

WRnE(*,*)’ ’ 

GOTO 20 
ENDIF 

IF ((rite.EQ.’Y’).OR.(rite.EQ.’y’)) THEN 
CALL Writefile(ktlpmod,bsIpmod,nurapix,lambda,RO,n) 

ENDIF 
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c DETERMINE THE AZIMUTHAL RMS PHASE ERROR 


30 WRITE(*,*)TIND THE AZIMUTHAL RMS PHASE ERROR? (Y/N)’ 
READ(*,*)mean 
WRITE(*,*)’ ’ 

IF ((mean.NE.’Y’)j\ND.(mean.NE.y)-AND. 

+ (mean.NE.’N’)j\ND.(mean.NE.’n’)) THEN 
WRITE(*,*)’ERROR, REENTER.’ 

WRITE(*,*)’ ’ 

GOTO 30 
ENDIF 

IF ((mean.EQ.’Y’).OR.(mean.EQ.’y’)) THEN 

c READ IN DATA 

CALL Readtrufile(TRUDATA,n) 

c DETERMINE SQUARE PHASE ERROR 

CALL AMSPE(KTDATA,BSDATA,TRUDATA, 

+ kterror,bseiTor,nyquist,n; 

c DETERMINE AZIMUTHAL AVERAGE OF THE RMS PHASE ERROR 

CALL AMSPEcalc(kterror,bserror,rKTerror,rBSeiTor, 

+ nyquist,n) 

c WRITE AMSPE TO A FILE FOR EACH CORRELATION TECHNIQUE 

CALL WriteiTfile(rKTeiTor,rBSerror,nyquist,lambda, 

+ RO,numpix,iktcount,itccount,n) 

ENDIF 

STOP 

END 
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SUBROUTINE LIST 


SUBROUTINE AMSPE(KTDAT^BSDATA,TRUDATA, 

+ kteiTor,bseiTor,nyqiiist,n) 

c THIS S/R CALCULATES THE SQUARE PHASE ERROR OF THE 

c RECONSTRUCTED PHASES 

COMPLEX* 16 KTDATA(n,n), BSDATA(n,n), TRUDATA(n,n) 
REAL*8 kterror(n,n), bseiTor(n,n), tniphase, ktphase, 

+ bsphase, pi, pi2 
pi = dacos(-1.0D+00) 
pi2 = 2 • pi 
n2pl = n/2 + 1 
DO 10 i = 1, n 
DO 10 j = 1, n 
X = float(j - (n2pl)) 
y = float((n2pl) - i) 
radius = sqrt(x**2.0 + y**2.0) 

IF (radius.L£.nyquist) THEN 
truphase = datan2(DIMAG(TRUDATA(ij)), 

+ DREAL(TRUDATA(ij))) 

ktphase = datan2(DIMAG(KTDATA(ij)), 

+ DREAL(KTDATA(y))) 

bsphase = datan2(DIMAG(BSDATA(ij)), 

+ DREAL(BSDATA(ij))) 

20 IF (truphase.GT.pi2) THEN 

truphase = truphase - pi2 
GOTO 20 
ENDIF 

30 IF (truphase.LT.-pi2) THEN 

truphase = truphase + pi2 
GOTO 30 
ENDIF 

40 IF (ktphase.GT.pi2) THEN 

ktphase = ktphase • pi2 
GOTO 40 
ENDIF 

50 IF (ktphase.LT-pi2) THEN 

ktphase = ktphase + pi2 
GOTO 50 
ENDIF 
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60 IF (bsphase.GT.pi2) THEN 
bsphase = bsphase - pi2 
GOTO 60 
ENDIF 

70 IF (bsphase.LT.-pi2) THEN 

bsphase = bsphase + pi2 
GOTO 70 
ENDIF 

kteiTor(ij) = truphase - ktphase 

80 IF (kterror(ij).GT.pi2) THEN 

kterror(ij) = kterror(ij) - pi2 
GOTO 80 
ENDIF 

90 IF (kterror(ij).LT.-pi2) THEN 

kterror(ij) = kterror(ij) + pi2 
GOTO 90 
ENDIF 

IF (kterror(ij).GT.pi) 

+ kteiTor(ij) = IrteiTor(ij) - pi2 

IF (kterror(ij).LT.-pi) 

+ kteiTor(ij) = Irterror(ij) + pi2 

kterror(ij) = (kteiTor(iJ))**2.0 

bseiTor(ij) = truphase - bsphase 

100 IF (bserror(iJ).GT.pi2) THEN 

bserror(ij) = bserror(ij) - pi2 
GOTO 100 
ENDIF 

no IF (bserror(ij).LT.-pi2) THEN 

bserror(ij) = bserror(ij) + pi2 
GOTO 110 
ENDIF 

IF (bserror(ij).GT.pi) 

+ bserror(ij) = bseiTor(ij) - pi2 

IF (bserror(iJ).LT.-pi) 

+ bserror(ij) = bseiTor(ij) + pi2 

bserror(ij) = (bserror(ij))**2.0 
ENDIF 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE AMSPEcalc(kterror,bserror,rKTerror, 

+ rBSerror^yquist,n) 

c THIS S/R CALCULATES THE RMS PHASE ERROR AS A FUNCTION OF 

c RADIUS 

REAL*8 kteiTor(n,n), bserror(n,n), rKTerror(ii/2), 

+ rBSerror(n/2) 

INTEGER r 
n2pl = 11 / 2+1 
DO 10 r = 0, i^quist 
nerr = 0 
DO 20 i = 1, n 
DO 20 j = 1, n 
X = float(j - (n2pl)) 
y = £loat((n2pl) - i) 
radius = sqrt(x**2.0 + y**2.0) 

IF ((radius.GE.float(r))AND. 

+ (radius.LT.float(r+l))) THEN 
nerr = nerr + 1 

rKTerror(r) = rKTerror(r) + kteiTor(ij) 
rBSerror(r) = rBSerror(r) + bserror(ij) 

ENDIF 

20 CONTINUE 

rKTerror(r) = dsqrt(rKTerror(r)/nerr) 
rBSerror(r) = dsqrt(rBSerror(r)^err) 

10 CONTINUE 
RETURN 
END 


SUBROUTINE Modulus(mod,DATA,n) 

c THIS S/R DETERMINES THE MODULUS OF A COMPLEX NUMBER 

COMPLEX*16 DATA(n,n) 

REAL*8 mod(n,n) 

DO 10 i = 1, n 
DO 10 j = 1, n 
mod(ij) = ABS(DATA(iJ)) 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE Normalize(data,n) 

c THIS S/R NORMALIZES THE OUTPUT DATA TO 1.0. 

REAL*8 data(n,n), maxval 
maxval = O.OD+00 
DO 10 i = 1, n 
DO 10 j = 1, n 

IF (data(iJ).GT.maxval) maxval = data(ij) 

10 CONTINUE 
DO 20 i = 1, n 
DO 20 j = 1, n 
data(ij) = data(ij)/maxval 
20 CONTINUE 
RETURN 
END 


SUBROUTINE Readfile(KTDATA,BSDATA,n) 
c THIS S/R READS THE RECONSTRUCTED DATA FROM A FILE 


10 


COMPLEX*16 KTDATA(n,n), BSDATA(n,n) 

CHARACTER* 16 ktfile, bsfile 
WRTTE(*,*)’ ’ 

WRrTE(*,*)’ENTER INPUT KT RECON FILE NAME (16 CHAR’ 
WRrTE(*,*)’MAX):’ 

READ (*,30) ktffle 
WRrTE(*,*)’ ’ 

WRITER,*)’ENTER INPUT TC RECON FILE NAME (16 CHAR’ 
WRrTE(*,*)’MAX):’ 

READ (*,30) bsffle 
WRITE(*,*)’ ’ 

OPEN(UNrT=40,FILE=ktfile,STATUS=’OLD’) 

OPEN(UNrT=50,FILE=bsfile,STATUS=’OLD’) 

DO 10 i = 1, n 


DO 10 j = 1, n 
READ(40,*) KTDATA(io) 
CONTINUE 
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DO 20 i = 1, n 
DO20j = l,n 
READ(50,*) BSDATA(ij) 
20 CONTINUE 
30 FORMAT(A16) 

RETURN 

END 


SUBROUTINE Readsiir(ktsnr,bssnr,ktr,bsr 4 iyquist,n) 

c THIS S/R READS THE SNR DATA FROM A FILE 

REAL*8 ktsnr(ii/2), bssnr(n/2) 

REAL ktr(n/2), bsr(n/2) 

CHARACTER* 16 ktsiirfile, bssnrfile 

WRITE(*,*)’ENTER KT SNR FILE NAME (16 CHAR MAX):’ 
READ (*.30) ktsnrfile 
WRnE(*,*)’ ’ 

WRITE(*,*)’ENTER TC SNR FILE NAME (16 CHAR MAX):’ 
READ (*,30) bssnrfile 
WRnE(*,*)’ ’ 

OPEN(UNIT=404TLE«ktsnrfile,STATUS=’OLD’) 

OPEN(UNIT»50,FILE«bssnrfilc,STATUS*’OLD’) 

DO 10 i = 1, nyquist 
READ(40,*) ktr(i), ktsnr(i) 

10 CONTINUE 

DO 20 i = 1, nyquist 
READ(50,*) bsr(i), bssnr(i) 

20 CONTINUE 
30 FORMAT(A16) 

RETURN 

END 


SUBROUTINE Readtrufile(TRUDATAn) 

c THIS S/R READS THE TRUTH DATA FROM A FILE < 

COMPLEX* 16 TRUDATA(n,n) 

CHARACTER* 16 file 

WRITE(*,*)’ENTER INPUT TRUTH DATA FILE NAME (16 CHAR’ 
WRITE(*,*)’MAX):’ 
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READ (*,20) file 
WRITE(*,*)’ ’ 

OPEN(UIOT=30JFILE=file,STATUS=’OLD’) 
DO 10 i = 1, n 
DO 10 j *= l,n 
READ(30,*) TRUDATA(ij) 

10 CONTINUE 
20 FORMAT(A16) 

RETURN 

END 


SUBROUTINE Truncate(DATA,LPDATA,trunc,lambda, 

+ R0,numpix4i) 

c THIS S/R FILTERS THE OBJECT SPECTRUM ARRAY BY 
c TRUNCATING THE ARRAY AT A RADIUS SET BY THE SNR VALUE 
c GREATER THAN 1 

COMPLEX* 16 DATA(n,n), LPDATA(n,n) 

REAL lambda 
n2 = 11/2 
n2pl = n2 + 1 
pi = acos(-1.0E+00) 

DO 10 i = 1, n 
DO 10 j = l,n 
X = float(j-n2pl) 
y = float(i-n2pl) 

rad = sqrt(x**2 + y**2) * numpix * (lambda/RO) * 

+ (180.0/pi) * 3600.0 

IF (rad.LE.trunc) THEN 
LPDATA(ij) = DATA(ij) 

ELSE 

LPDATA(ij) = (0.0,0.0) 

ENDIF 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE Truncval(snr,r,trunc,nyquist,icount,n) 


c THIS S/R DETERMINES THE TRUNCATION VALUE BASED ON THE 
c SNR. THE VALUE IS BASED ON THE POINT AT WHICH THE SNR IS 
c UNITY OR GREATER. 

REAL*8 snr(n/2) 

REALr(n/2) 

j = 1 

DO 10 i = 1, nyquist 

IF ((snr(i).GE.1.0D+00).AND.O.EQ.i)) THEN 
trunc = r(i) 
icount = j 
j = j + 1 
ENDIF 

10 CONTINUE 
RETURN 
END 


SUBROUTINE Writefile(ktlpdata,bslpdata,numpix, 
+ lambda,R0,n) 

c THIS S/R WRITES THE DATA TO A FILE 


REAL*8 ktlpdata(n,n), bslpdata(n,n) 

REAL conv, lambda, RO, x, y 
CHARACTER* 16 ktlpfile, bslpflle 
pi = acos(-1.0E+00) 

WRITE(*,*)’ENTER OUTPUT KT FILE NAME (16 CHAR MAX):’ 
READ (*,50) ktlpfile 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER OUTPUT TC FILE NAME (16 CHAR MAX):’ 
READ (*,50) bslpfile 
WRITE(*,*)’ ’ 

OPEN(UNIT=30,FILE=ktlpfile,STATUS=’NEW’) 
OPEN(UNIT=40,FILE=bslpfile,STATUS=’NEW’) 

DO 10 i = 1, n 
DO 10 j - 1, n 

conv = float(numpix) * (lambda/RO) * 

+ (180.0/pi) * 3600.0 

X = float(j - (n/z+l)) * conv/float(n) 
y = float(i - (n/2+1)) * conv/float(n) 
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WRITE(30,60) X, y, ktlpdata(y) 

10 CX)NTINIJE 
DO 20 i = 1, n 
DO 20 j = 1, n 

conv = float(numpix) * (lambda/RO) * 
+ (180.0/pi) • 3600.0 

X = float(} - (n^+1)) • conv/float(n) 
y = float(i - (n/2+1)) • conv/float(n) 
WRITE(40,60) X, y, bslpdata(ij) 

20 CONTINUE 
50 FORMAT(A16) 

60 FORMAT(F7.4,2x,F7.4,2x,F7.4) 

RETURN 
END 


SUBROUTINE Writerrfile(rKTerror,rBSeiTor,nyquist, 

+ lambda,RO,numpix,iktcount, 

+ itccount,n) 

c THIS S/R WRITES THE ERROR DATA TO A FILE 

REAL*8 rKTeiTor(ii/2), rBSeiTor(nA2) 

REAL X, y, lambda 
CHARACTER* 16 ktfile, bsfile 
INTEGER r 
pi = acos(-1.0+00) 

WRITE(*,*)’ENTER OUTPUT KT ERROR FILE NAME (16 CHAR’ 
WRITE(*,*)’MAX);’ 

READ (*,30) ktfile 
WRITE(*,*)’ ’ 

WRITE(*,*)’ENTER OUTPUT TC ERROR FILE NAME (16 CHAR’ 
WRITE(*,*)’MAX):’ 

READ (*,30) bsfile 
WRITE(*,*)’ ’ 

OPEN(UNIT= l,FILE=ktfile,STATUS=’NEW’) 
OPEN(UNIT=2,FILE=bsfile,STATUS=’NEW’) 

DO 10 r = 1, iktcount 

X = r * numpix * (lambda/RO) * (180.0/pi) * 3600.0 
WRrrE(l,*) X, rKTerror(r) 

10 CONTINUE 
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DO 20 r = 1, itccount 

X = r * numpix * (lambda/RO) * (180.0/pi) * 3600.0 
WRITE(2,*) X, rBSerror(r) 

20 CONTINUE 
30 FORMAT(A16) 

RETURN 

END 
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